src/INTERP_KERNEL/VolSurfUser.txx

   1 // Copyright (C) 2007-2024  CEA, EDF
   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 #ifndef __VOLSURFUSER_TXX__
  21 #define __VOLSURFUSER_TXX__
  22
  23 #include "VolSurfUser.hxx"
  24 #include "VolSurfFormulae.hxx"
  25 #include "InterpolationUtils.hxx"
  26 #include "VectorUtils.hxx"
  27
  28 #include <algorithm>
  29
  30 namespace INTERP_KERNEL
  31 {
  32   template<class ConnType, NumberingPolicy numPol, int SPACEDIM>
  33   double computeVolSurfOfCell(NormalizedCellType type, const ConnType *connec, mcIdType lgth, const double *coords)
  34   {
  35     switch(type)
  36       {
  37       case INTERP_KERNEL::NORM_SEG2 :
  38       case INTERP_KERNEL::NORM_SEG4 :
  39         {
  40           ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
  41           ConnType 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_SEG3 :
  45         {
  46           ConnType beginNode = OTT<ConnType,numPol>::coo2C(connec[0]);
  47           ConnType endNode = OTT<ConnType,numPol>::coo2C(connec[1]);
  48           ConnType middleNode = OTT<ConnType,numPol>::coo2C(connec[2]);
  49           return INTERP_KERNEL::calculateLgthForSeg3(coords+(SPACEDIM*beginNode),coords+(SPACEDIM*endNode),coords+(SPACEDIM*middleNode),SPACEDIM);
  50         }
  51       case INTERP_KERNEL::NORM_TRI3 :
  52         {
  53           ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
  54           ConnType N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
  55           ConnType N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
  56
  57           return INTERP_KERNEL::calculateAreaForTria(coords+(SPACEDIM*N1),
  58                                                      coords+(SPACEDIM*N2),
  59                                                      coords+(SPACEDIM*N3),
  60                                                      SPACEDIM);
  61         }
  62         break;
  63
  64       case INTERP_KERNEL::NORM_TRI6 :
  65       case INTERP_KERNEL::NORM_TRI7 :
  66         {
  67           const double *pts[6];
  68           pts[0] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]);
  69           pts[1] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]);
  70           pts[2] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]);
  71           pts[3] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[3]);
  72           pts[4] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[4]);
  73           pts[5] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[5]);
  74           return INTERP_KERNEL::calculateAreaForQPolyg(pts,6,SPACEDIM);
  75         }
  76         break;
  77       case INTERP_KERNEL::NORM_QUAD4 :
  78         {
  79           ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
  80           ConnType N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
  81           ConnType N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
  82           ConnType N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
  83
  84           return INTERP_KERNEL::calculateAreaForQuad(coords+SPACEDIM*N1,
  85                                                      coords+SPACEDIM*N2,
  86                                                      coords+SPACEDIM*N3,
  87                                                      coords+SPACEDIM*N4,
  88                                                      SPACEDIM);
  89         }
  90         break;
  91       case INTERP_KERNEL::NORM_QUAD8 :
  92       case INTERP_KERNEL::NORM_QUAD9 :
  93         {
  94           const double *pts[8];
  95           pts[0] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]);
  96           pts[1] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]);
  97           pts[2] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]);
  98           pts[3] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[3]);
  99           pts[4] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[4]);
 100           pts[5] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[5]);
 101           pts[6] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[6]);
 102           pts[7] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[7]);
 103           return INTERP_KERNEL::calculateAreaForQPolyg(pts,8,SPACEDIM);
 104         }
 105         break;
 106       case INTERP_KERNEL::NORM_POLYGON :
 107         {
 108           const double **pts=new const double *[lgth];
 109           for(int inod=0;inod<lgth;inod++)
 110             pts[inod] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[inod]);
 111           double val=INTERP_KERNEL::calculateAreaForPolyg(pts,lgth,SPACEDIM);
 112           delete [] pts;
 113           return val;
 114         }
 115         break;
 116       case INTERP_KERNEL::NORM_QPOLYG :
 117         {
 118           const double **pts=new const double *[lgth];
 119           for(int inod=0;inod<lgth;inod++)
 120             pts[inod] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[inod]);
 121           double val=INTERP_KERNEL::calculateAreaForQPolyg(pts,lgth,SPACEDIM);
 122           delete [] pts;
 123           return val;
 124         }
 125         break;
 126       case INTERP_KERNEL::NORM_TETRA4 :
 127       case INTERP_KERNEL::NORM_TETRA10 :
 128         {
 129           ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
 130           ConnType N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
 131           ConnType N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
 132           ConnType N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
 133
 134           return INTERP_KERNEL::calculateVolumeForTetra(coords+SPACEDIM*N1,
 135                                                         coords+SPACEDIM*N2,
 136                                                         coords+SPACEDIM*N3,
 137                                                         coords+SPACEDIM*N4);
 138         }
 139         break;
 140
 141       case INTERP_KERNEL::NORM_PYRA5 :
 142       case INTERP_KERNEL::NORM_PYRA13 :
 143         {
 144           ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
 145           ConnType N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
 146           ConnType N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
 147           ConnType N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
 148           ConnType N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
 149
 150           return INTERP_KERNEL::calculateVolumeForPyra(coords+SPACEDIM*N1,
 151                                                        coords+SPACEDIM*N2,
 152                                                        coords+SPACEDIM*N3,
 153                                                        coords+SPACEDIM*N4,
 154                                                        coords+SPACEDIM*N5);
 155         }
 156         break;
 157
 158       case INTERP_KERNEL::NORM_PENTA6 :
 159       case INTERP_KERNEL::NORM_PENTA15 :
 160       case INTERP_KERNEL::NORM_PENTA18 :
 161         {
 162           ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
 163           ConnType N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
 164           ConnType N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
 165           ConnType N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
 166           ConnType N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
 167           ConnType N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
 168
 169           return INTERP_KERNEL::calculateVolumeForPenta(coords+SPACEDIM*N1,
 170                                                         coords+SPACEDIM*N2,
 171                                                         coords+SPACEDIM*N3,
 172                                                         coords+SPACEDIM*N4,
 173                                                         coords+SPACEDIM*N5,
 174                                                         coords+SPACEDIM*N6);
 175         }
 176         break;
 177
 178       case INTERP_KERNEL::NORM_HEXA8 :
 179       case INTERP_KERNEL::NORM_HEXA20 :
 180       case INTERP_KERNEL::NORM_HEXA27 :
 181         {
 182           ConnType N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
 183           ConnType N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
 184           ConnType N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
 185           ConnType N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
 186           ConnType N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
 187           ConnType N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
 188           ConnType N7 = OTT<ConnType,numPol>::coo2C(connec[6]);
 189           ConnType N8 = OTT<ConnType,numPol>::coo2C(connec[7]);
 190
 191           return INTERP_KERNEL::calculateVolumeForHexa(coords+SPACEDIM*N1,
 192                                                        coords+SPACEDIM*N2,
 193                                                        coords+SPACEDIM*N3,
 194                                                        coords+SPACEDIM*N4,
 195                                                        coords+SPACEDIM*N5,
 196                                                        coords+SPACEDIM*N6,
 197                                                        coords+SPACEDIM*N7,
 198                                                        coords+SPACEDIM*N8);
 199         }
 200         break;
 201       case INTERP_KERNEL::NORM_HEXGP12:
 202         {
 203           const ConnType connecHexa12[43]={
 204             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,
 205             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,
 206             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,
 207             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,
 208             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,
 209             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,
 210             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,
 211             OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[0])};
 212           return calculateVolumeForPolyh2<ConnType,numPol>(connecHexa12,43,coords);
 213         }
 214       case INTERP_KERNEL::NORM_POLYHED :
 215         {
 216           return calculateVolumeForPolyh2<ConnType,numPol>(connec,lgth,coords);
 217         }
 218         break;
 219       default:
 220         throw INTERP_KERNEL::Exception("Not recognized cell type to get Length/Area/Volume on it !");
 221       }
 222   }
 223
 224   template<class ConnType, NumberingPolicy numPolConn>
 225   double computeVolSurfOfCell2(NormalizedCellType type, const ConnType *connec, mcIdType lgth, const double *coords, int spaceDim)
 226   {
 227     if(spaceDim==3)
 228       return computeVolSurfOfCell<ConnType,numPolConn,3>(type,connec,lgth,coords);
 229     if(spaceDim==2)
 230       return computeVolSurfOfCell<ConnType,numPolConn,2>(type,connec,lgth,coords);
 231     if(spaceDim==1)
 232       return computeVolSurfOfCell<ConnType,numPolConn,1>(type,connec,lgth,coords);
 233     throw INTERP_KERNEL::Exception("Invalid spaceDim specified : must be 1, 2 or 3");
 234   }
 235
 236
 237   template<class ConnType, NumberingPolicy numPol,int SPACEDIM>
 238   void computeBarycenter(NormalizedCellType type, const ConnType *connec, mcIdType lgth, const double *coords, double *res)
 239   {
 240     switch(type)
 241       {
 242       case NORM_SEG2:
 243       case NORM_SEG4:
 244         {
 245           std::copy(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]),
 246                     coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]+1),res);
 247           std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]),res,std::plus<double>());
 248           std::transform(res,res+SPACEDIM,res,std::bind(std::multiplies<double>(),std::placeholders::_1,0.5));
 249           break;
 250         }
 251       case NORM_SEG3:
 252         {
 253           if(SPACEDIM==2)
 254             {
 255               Edge *ed=Edge::BuildEdgeFrom3Points(coords+2*OTT<ConnType,numPol>::coo2C(connec[0]),coords+2*OTT<ConnType,numPol>::coo2C(connec[2]),coords+2*OTT<ConnType,numPol>::coo2C(connec[1]));
 256               ed->getBarycenter(res);
 257               ed->decrRef();
 258             }
 259           else if(SPACEDIM==1)
 260             {
 261               *res=(coords[OTT<ConnType,numPol>::coo2C(connec[0])]+coords[OTT<ConnType,numPol>::coo2C(connec[1])])/2.;
 262             }
 263           else if(SPACEDIM==3)
 264             {
 265               std::copy(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]),
 266                         coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]+1),res);
 267               std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]),res,std::plus<double>());
 268               std::transform(res,res+SPACEDIM,res,std::bind(std::multiplies<double>(),std::placeholders::_1,0.5));
 269             }
 270           else
 271             throw INTERP_KERNEL::Exception("computeBarycenter for SEG3 only SPACEDIM 1,2 or 3 supported !");
 272           break;
 273         }
 274       case NORM_TRI3:
 275       case NORM_TRI7:
 276         {
 277           std::copy(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]),
 278                     coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]+1),res);
 279           std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]),res,std::plus<double>());
 280           std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]),res,std::plus<double>());
 281           std::transform(res,res+SPACEDIM,res,std::bind(std::multiplies<double>(),std::placeholders::_1,1./3.));
 282           break;
 283         }
 284       case NORM_TRI6:
 285         {
 286           if(SPACEDIM==2)
 287             {
 288               double *pts[6];
 289               pts[0] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]));
 290               pts[1] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]));
 291               pts[2] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]));
 292               pts[3] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[3]));
 293               pts[4] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[4]));
 294               pts[5] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[5]));
 295               computeQPolygonBarycenter2D(pts,6,2,res);
 296             }
 297           else if(SPACEDIM==3)
 298             computePolygonBarycenter3D<ConnType,numPol>(connec,lgth/2,coords,res);
 299           else
 300             throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
 301           break;
 302         }
 303       case NORM_QUAD4:
 304       case NORM_POLYGON:
 305         {
 306           if(SPACEDIM==2)
 307             computePolygonBarycenter2D<ConnType,numPol>(connec,lgth,coords,res);
 308           else if(SPACEDIM==3)
 309             computePolygonBarycenter3D<ConnType,numPol>(connec,lgth,coords,res);
 310           else
 311             throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
 312           break;
 313         }
 314       case NORM_QUAD8:
 315         {
 316           if(SPACEDIM==2)
 317             {
 318               double *pts[8];
 319               pts[0] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]));
 320               pts[1] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]));
 321               pts[2] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]));
 322               pts[3] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[3]));
 323               pts[4] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[4]));
 324               pts[5] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[5]));
 325               pts[6] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[6]));
 326               pts[7] = const_cast<double *>(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[7]));
 327               computeQPolygonBarycenter2D(pts,8,2,res);
 328             }
 329           else if(SPACEDIM==3)
 330             computePolygonBarycenter3D<ConnType,numPol>(connec,lgth/2,coords,res);
 331           else
 332             throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
 333           break;
 334         }
 335       case INTERP_KERNEL::NORM_QPOLYG :
 336         {
 337           if(SPACEDIM==2)
 338             {
 339               double **pts=new double *[lgth];
 340               for(int i=0;i<lgth;i++)
 341                 pts[i]=const_cast<double *>(coords+2*OTT<ConnType,numPol>::coo2C(connec[i]));
 342               computeQPolygonBarycenter2D(pts,lgth,2,res);
 343               delete [] pts;
 344             }
 345           else if(SPACEDIM==3)
 346             computePolygonBarycenter3D<ConnType,numPol>(connec,lgth/2,coords,res);
 347           else
 348             throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
 349           break;
 350         }
 351         break;
 352       case NORM_TETRA4:
 353         {
 354           res[0]=coords[3*OTT<ConnType,numPol>::coo2C(connec[0])];
 355           res[1]=coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+1];
 356           res[2]=coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+2];
 357           res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[1])];
 358           res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[1])+1];
 359           res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[1])+2];
 360           res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[2])];
 361           res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[2])+1];
 362           res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[2])+2];
 363           res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[3])];
 364           res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[3])+1];
 365           res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[3])+2];
 366           res[0]*=0.25; res[1]*=0.25; res[2]*=0.25;
 367           break;
 368         }
 369       case NORM_PYRA5:
 370         {
 371           double tmp[3];
 372           computePolygonBarycenter3D<ConnType,numPol>(connec,lgth-1,coords,tmp);
 373           res[0]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[4])]+3.*tmp[0])/4.;
 374           res[1]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[4])+1]+3.*tmp[1])/4.;
 375           res[2]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[4])+2]+3.*tmp[2])/4.;
 376           break;
 377         }
 378       case NORM_HEXA8:
 379         {
 380           const ConnType conn[29]={
 381             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,
 382             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,
 383             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,
 384             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,
 385             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,
 386             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[4]),OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[1]),
 387             };
 388           barycenterOfPolyhedron<ConnType,numPol>(conn,29,coords,res);
 389           break;
 390         }
 391       case NORM_PENTA6:
 392         {
 393           const ConnType conn[22]={
 394             OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[2]),-1,
 395             OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[4]),-1,
 396             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,
 397             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,
 398             OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[4])
 399           };
 400           barycenterOfPolyhedron<ConnType,numPol>(conn,22,coords,res);
 401           break;
 402         }
 403       case INTERP_KERNEL::NORM_HEXGP12:
 404         {
 405           const ConnType connecHexa12[43]={
 406             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,
 407             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,
 408             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,
 409             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,
 410             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,
 411             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,
 412             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,
 413             OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[0])};
 414           barycenterOfPolyhedron<ConnType,numPol>(connecHexa12,43,coords,res);
 415           break;
 416         }
 417       case NORM_POLYHED:
 418         {
 419           barycenterOfPolyhedron<ConnType,numPol>(connec,lgth,coords,res);
 420           break;
 421         }
 422       default:
 423         throw INTERP_KERNEL::Exception("Not recognized cell type to get Barycenter on it !");
 424       }
 425   }
 426
 427   template<class ConnType, NumberingPolicy numPolConn>
 428   void computeBarycenter2(NormalizedCellType type, const ConnType *connec, mcIdType lgth, const double *coords, int spaceDim, double *res)
 429   {
 430     if(spaceDim==3)
 431       return computeBarycenter<ConnType,numPolConn,3>(type,connec,lgth,coords,res);
 432     if(spaceDim==2)
 433       return computeBarycenter<ConnType,numPolConn,2>(type,connec,lgth,coords,res);
 434     if(spaceDim==1)
 435       return computeBarycenter<ConnType,numPolConn,1>(type,connec,lgth,coords,res);
 436     throw INTERP_KERNEL::Exception("Invalid spaceDim specified for compute barycenter : must be 1, 2 or 3");
 437   }
 438
 439   template<int SPACEDIM>
 440   void ComputeTriangleHeight(const double *PA, const double *PB, const double *PC, double *res)
 441   {
 442     constexpr double EPS = 1e-12;
 443     double AB = getDistanceBtw2Pts<SPACEDIM>(PA,PB); double maxLength = AB;
 444     double BC = getDistanceBtw2Pts<SPACEDIM>(PB,PC); if (BC > maxLength) maxLength = BC;
 445     double CA = getDistanceBtw2Pts<SPACEDIM>(PC,PA); if (CA > maxLength) maxLength = CA;
 446     double perim( (AB+BC+CA)*0.5 );
 447     double num( 2*sqrt(perim*(perim-AB)*(perim-BC)*(perim-CA)) );
 448     res[0] = num/AB; res[1] = num/BC; res[2] = num/CA;
 449     if (AB/maxLength <= EPS)
 450       res[0] = BC;
 451     if (BC/maxLength <= EPS)
 452       res[1] = CA;
 453     if (CA/maxLength <= EPS)
 454       res[3] = AB;
 455   }
 456 }
 457
 458 #endif