1 // Copyright (C) 2007-2012 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 :
45 case INTERP_KERNEL::NORM_TRI6 :
46 case INTERP_KERNEL::NORM_TRI7 :
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]);
52 return INTERP_KERNEL::calculateAreaForTria(coords+(SPACEDIM*N1),
59 case INTERP_KERNEL::NORM_QUAD4 :
60 case INTERP_KERNEL::NORM_QUAD8 :
61 case INTERP_KERNEL::NORM_QUAD9 :
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]);
68 return INTERP_KERNEL::calculateAreaForQuad(coords+SPACEDIM*N1,
76 case INTERP_KERNEL::NORM_POLYGON :
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);
86 case INTERP_KERNEL::NORM_TETRA4 :
87 case INTERP_KERNEL::NORM_TETRA10 :
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]);
94 return INTERP_KERNEL::calculateVolumeForTetra(coords+SPACEDIM*N1,
101 case INTERP_KERNEL::NORM_PYRA5 :
102 case INTERP_KERNEL::NORM_PYRA13 :
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]);
110 return INTERP_KERNEL::calculateVolumeForPyra(coords+SPACEDIM*N1,
118 case INTERP_KERNEL::NORM_PENTA6 :
119 case INTERP_KERNEL::NORM_PENTA15 :
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]);
128 return INTERP_KERNEL::calculateVolumeForPenta(coords+SPACEDIM*N1,
137 case INTERP_KERNEL::NORM_HEXA8 :
138 case INTERP_KERNEL::NORM_HEXA20 :
139 case INTERP_KERNEL::NORM_HEXA27 :
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]);
150 return INTERP_KERNEL::calculateVolumeForHexa(coords+SPACEDIM*N1,
160 case INTERP_KERNEL::NORM_HEXGP12:
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);
173 case INTERP_KERNEL::NORM_POLYHED :
175 return calculateVolumeForPolyh2<ConnType,numPol>(connec,lgth,coords);
179 throw INTERP_KERNEL::Exception("Not recognized cell type to get Length/Area/Volume on it !");
183 template<class ConnType, NumberingPolicy numPolConn>
184 double computeVolSurfOfCell2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim)
187 return computeVolSurfOfCell<ConnType,numPolConn,3>(type,connec,lgth,coords);
189 return computeVolSurfOfCell<ConnType,numPolConn,2>(type,connec,lgth,coords);
191 return computeVolSurfOfCell<ConnType,numPolConn,1>(type,connec,lgth,coords);
192 throw INTERP_KERNEL::Exception("Invalid spaceDim specified : must be 1, 2 or 3");
196 template<class ConnType, NumberingPolicy numPol,int SPACEDIM>
197 void computeBarycenter(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, double *res)
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));
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.));
226 computePolygonBarycenter2D<ConnType,numPol>(connec,lgth,coords,res);
228 computePolygonBarycenter3D<ConnType,numPol>(connec,lgth,coords,res);
230 throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
236 computePolygonBarycenter2D<ConnType,numPol>(connec,lgth/2,coords,res);
238 computePolygonBarycenter3D<ConnType,numPol>(connec,lgth/2,coords,res);
240 throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
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;
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.;
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]),
279 barycenterOfPolyhedron<ConnType,numPol>(conn,29,coords,res);
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])
291 barycenterOfPolyhedron<ConnType,numPol>(conn,22,coords,res);
294 case INTERP_KERNEL::NORM_HEXGP12:
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);
310 barycenterOfPolyhedron<ConnType,numPol>(connec,lgth,coords,res);
314 throw INTERP_KERNEL::Exception("Not recognized cell type to get Barycenter on it !");
318 template<class ConnType, NumberingPolicy numPolConn>
319 void computeBarycenter2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim, double *res)
322 return computeBarycenter<ConnType,numPolConn,3>(type,connec,lgth,coords,res);
324 return computeBarycenter<ConnType,numPolConn,2>(type,connec,lgth,coords,res);
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");