1 // Copyright (C) 2007-2013 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
20 #ifndef __VOLSURFUSER_TXX__
21 #define __VOLSURFUSER_TXX__
23 #include "VolSurfUser.hxx"
24 #include "VolSurfFormulae.hxx"
25 #include "InterpolationUtils.hxx"
29 namespace INTERP_KERNEL
31 template<class ConnType, NumberingPolicy numPol, int SPACEDIM>
32 double computeVolSurfOfCell(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords)
36 case INTERP_KERNEL::NORM_SEG2 :
37 case INTERP_KERNEL::NORM_SEG3 :
38 case INTERP_KERNEL::NORM_SEG4 :
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);
44 case INTERP_KERNEL::NORM_TRI3 :
46 int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
47 int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
48 int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
50 return INTERP_KERNEL::calculateAreaForTria(coords+(SPACEDIM*N1),
57 case INTERP_KERNEL::NORM_TRI6 :
58 case INTERP_KERNEL::NORM_TRI7 :
61 pts[0] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]);
62 pts[1] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]);
63 pts[2] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]);
64 pts[3] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[3]);
65 pts[4] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[4]);
66 pts[5] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[5]);
67 return INTERP_KERNEL::calculateAreaForQPolyg(pts,6,SPACEDIM);
70 case INTERP_KERNEL::NORM_QUAD4 :
72 int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
73 int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
74 int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
75 int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
77 return INTERP_KERNEL::calculateAreaForQuad(coords+SPACEDIM*N1,
84 case INTERP_KERNEL::NORM_QUAD8 :
85 case INTERP_KERNEL::NORM_QUAD9 :
88 pts[0] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]);
89 pts[1] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]);
90 pts[2] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]);
91 pts[3] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[3]);
92 pts[4] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[4]);
93 pts[5] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[5]);
94 pts[6] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[6]);
95 pts[7] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[7]);
96 return INTERP_KERNEL::calculateAreaForQPolyg(pts,8,SPACEDIM);
99 case INTERP_KERNEL::NORM_POLYGON :
101 const double **pts=new const double *[lgth];
102 for(int inod=0;inod<lgth;inod++)
103 pts[inod] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[inod]);
104 double val=INTERP_KERNEL::calculateAreaForPolyg(pts,lgth,SPACEDIM);
109 case INTERP_KERNEL::NORM_QPOLYG :
111 const double **pts=new const double *[lgth];
112 for(int inod=0;inod<lgth;inod++)
113 pts[inod] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[inod]);
114 double val=INTERP_KERNEL::calculateAreaForQPolyg(pts,lgth,SPACEDIM);
118 case INTERP_KERNEL::NORM_TETRA4 :
119 case INTERP_KERNEL::NORM_TETRA10 :
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]);
126 return INTERP_KERNEL::calculateVolumeForTetra(coords+SPACEDIM*N1,
133 case INTERP_KERNEL::NORM_PYRA5 :
134 case INTERP_KERNEL::NORM_PYRA13 :
136 int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
137 int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
138 int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
139 int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
140 int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
142 return INTERP_KERNEL::calculateVolumeForPyra(coords+SPACEDIM*N1,
150 case INTERP_KERNEL::NORM_PENTA6 :
151 case INTERP_KERNEL::NORM_PENTA15 :
153 int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
154 int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
155 int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
156 int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
157 int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
158 int N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
160 return INTERP_KERNEL::calculateVolumeForPenta(coords+SPACEDIM*N1,
169 case INTERP_KERNEL::NORM_HEXA8 :
170 case INTERP_KERNEL::NORM_HEXA20 :
171 case INTERP_KERNEL::NORM_HEXA27 :
173 int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
174 int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
175 int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
176 int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
177 int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
178 int N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
179 int N7 = OTT<ConnType,numPol>::coo2C(connec[6]);
180 int N8 = OTT<ConnType,numPol>::coo2C(connec[7]);
182 return INTERP_KERNEL::calculateVolumeForHexa(coords+SPACEDIM*N1,
192 case INTERP_KERNEL::NORM_HEXGP12:
194 const int connecHexa12[43]={
195 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,
196 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,
197 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,
198 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,
199 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,
200 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,
201 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,
202 OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[0])};
203 return calculateVolumeForPolyh2<ConnType,numPol>(connecHexa12,43,coords);
205 case INTERP_KERNEL::NORM_POLYHED :
207 return calculateVolumeForPolyh2<ConnType,numPol>(connec,lgth,coords);
211 throw INTERP_KERNEL::Exception("Not recognized cell type to get Length/Area/Volume on it !");
215 template<class ConnType, NumberingPolicy numPolConn>
216 double computeVolSurfOfCell2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim)
219 return computeVolSurfOfCell<ConnType,numPolConn,3>(type,connec,lgth,coords);
221 return computeVolSurfOfCell<ConnType,numPolConn,2>(type,connec,lgth,coords);
223 return computeVolSurfOfCell<ConnType,numPolConn,1>(type,connec,lgth,coords);
224 throw INTERP_KERNEL::Exception("Invalid spaceDim specified : must be 1, 2 or 3");
228 template<class ConnType, NumberingPolicy numPol,int SPACEDIM>
229 void computeBarycenter(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, double *res)
237 std::copy(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]),
238 coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]+1),res);
239 std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]),res,std::plus<double>());
240 std::transform(res,res+SPACEDIM,res,std::bind2nd(std::multiplies<double>(),0.5));
247 std::copy(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]),
248 coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]+1),res);
249 std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]),res,std::plus<double>());
250 std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]),res,std::plus<double>());
251 std::transform(res,res+SPACEDIM,res,std::bind2nd(std::multiplies<double>(),1./3.));
258 computePolygonBarycenter2D<ConnType,numPol>(connec,lgth,coords,res);
260 computePolygonBarycenter3D<ConnType,numPol>(connec,lgth,coords,res);
262 throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
268 computePolygonBarycenter2D<ConnType,numPol>(connec,lgth/2,coords,res);
270 computePolygonBarycenter3D<ConnType,numPol>(connec,lgth/2,coords,res);
272 throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
277 res[0]=coords[3*OTT<ConnType,numPol>::coo2C(connec[0])];
278 res[1]=coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+1];
279 res[2]=coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+2];
280 res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[1])];
281 res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[1])+1];
282 res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[1])+2];
283 res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[2])];
284 res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[2])+1];
285 res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[2])+2];
286 res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[3])];
287 res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[3])+1];
288 res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[3])+2];
289 res[0]*=0.25; res[1]*=0.25; res[2]*=0.25;
295 computePolygonBarycenter3D<ConnType,numPol>(connec,lgth-1,coords,tmp);
296 res[0]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[4])]+3.*tmp[0])/4.;
297 res[1]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[4])+1]+3.*tmp[1])/4.;
298 res[2]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[4])+2]+3.*tmp[2])/4.;
304 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,
305 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,
306 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,
307 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,
308 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,
309 OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[4]),OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[1]),
311 barycenterOfPolyhedron<ConnType,numPol>(conn,29,coords,res);
317 OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[2]),-1,
318 OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[4]),-1,
319 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,
320 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,
321 OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[4])
323 barycenterOfPolyhedron<ConnType,numPol>(conn,22,coords,res);
326 case INTERP_KERNEL::NORM_HEXGP12:
328 const int connecHexa12[43]={
329 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,
330 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,
331 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,
332 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,
333 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,
334 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,
335 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,
336 OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[0])};
337 barycenterOfPolyhedron<ConnType,numPol>(connecHexa12,43,coords,res);
342 barycenterOfPolyhedron<ConnType,numPol>(connec,lgth,coords,res);
346 throw INTERP_KERNEL::Exception("Not recognized cell type to get Barycenter on it !");
350 template<class ConnType, NumberingPolicy numPolConn>
351 void computeBarycenter2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim, double *res)
354 return computeBarycenter<ConnType,numPolConn,3>(type,connec,lgth,coords,res);
356 return computeBarycenter<ConnType,numPolConn,2>(type,connec,lgth,coords,res);
358 return computeBarycenter<ConnType,numPolConn,1>(type,connec,lgth,coords,res);
359 throw INTERP_KERNEL::Exception("Invalid spaceDim specified for compute barycenter : must be 1, 2 or 3");