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