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 #ifndef __VOLSURFUSER_TXX__
20 #define __VOLSURFUSER_TXX__
22 #include "VolSurfUser.hxx"
23 #include "VolSurfFormulae.hxx"
24 #include "InterpolationUtils.hxx"
28 namespace INTERP_KERNEL
30 template<class ConnType, NumberingPolicy numPol, int SPACEDIM>
31 double computeVolSurfOfCell(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords)
35 case INTERP_KERNEL::NORM_SEG2 :
36 case INTERP_KERNEL::NORM_SEG3 :
37 case INTERP_KERNEL::NORM_SEG4 :
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);
43 case INTERP_KERNEL::NORM_TRI3 :
44 case INTERP_KERNEL::NORM_TRI6 :
45 case INTERP_KERNEL::NORM_TRI7 :
47 int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
48 int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
49 int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
51 return INTERP_KERNEL::calculateAreaForTria(coords+(SPACEDIM*N1),
58 case INTERP_KERNEL::NORM_QUAD4 :
59 case INTERP_KERNEL::NORM_QUAD8 :
60 case INTERP_KERNEL::NORM_QUAD9 :
62 int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
63 int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
64 int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
65 int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
67 return INTERP_KERNEL::calculateAreaForQuad(coords+SPACEDIM*N1,
75 case INTERP_KERNEL::NORM_POLYGON :
77 const double **pts=new const double *[lgth];
78 for(int inod=0;inod<lgth;inod++)
79 pts[inod] = coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[inod]);
80 double val=INTERP_KERNEL::calculateAreaForPolyg(pts,lgth,SPACEDIM);
85 case INTERP_KERNEL::NORM_TETRA4 :
86 case INTERP_KERNEL::NORM_TETRA10 :
88 int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
89 int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
90 int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
91 int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
93 return INTERP_KERNEL::calculateVolumeForTetra(coords+SPACEDIM*N1,
100 case INTERP_KERNEL::NORM_PYRA5 :
101 case INTERP_KERNEL::NORM_PYRA13 :
103 int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
104 int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
105 int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
106 int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
107 int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
109 return INTERP_KERNEL::calculateVolumeForPyra(coords+SPACEDIM*N1,
117 case INTERP_KERNEL::NORM_PENTA6 :
118 case INTERP_KERNEL::NORM_PENTA15 :
120 int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
121 int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
122 int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
123 int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
124 int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
125 int N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
127 return INTERP_KERNEL::calculateVolumeForPenta(coords+SPACEDIM*N1,
136 case INTERP_KERNEL::NORM_HEXA8 :
137 case INTERP_KERNEL::NORM_HEXA20 :
138 case INTERP_KERNEL::NORM_HEXA27 :
140 int N1 = OTT<ConnType,numPol>::coo2C(connec[0]);
141 int N2 = OTT<ConnType,numPol>::coo2C(connec[1]);
142 int N3 = OTT<ConnType,numPol>::coo2C(connec[2]);
143 int N4 = OTT<ConnType,numPol>::coo2C(connec[3]);
144 int N5 = OTT<ConnType,numPol>::coo2C(connec[4]);
145 int N6 = OTT<ConnType,numPol>::coo2C(connec[5]);
146 int N7 = OTT<ConnType,numPol>::coo2C(connec[6]);
147 int N8 = OTT<ConnType,numPol>::coo2C(connec[7]);
149 return INTERP_KERNEL::calculateVolumeForHexa(coords+SPACEDIM*N1,
159 case INTERP_KERNEL::NORM_HEXGP12:
161 const int connecHexa12[43]={
162 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,
163 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,
164 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,
165 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,
166 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,
167 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,
168 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,
169 OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[0])};
170 return calculateVolumeForPolyh2<ConnType,numPol>(connecHexa12,43,coords);
172 case INTERP_KERNEL::NORM_POLYHED :
174 return calculateVolumeForPolyh2<ConnType,numPol>(connec,lgth,coords);
178 throw INTERP_KERNEL::Exception("Not recognized cell type to get Length/Area/Volume on it !");
182 template<class ConnType, NumberingPolicy numPolConn>
183 double computeVolSurfOfCell2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim)
186 return computeVolSurfOfCell<ConnType,numPolConn,3>(type,connec,lgth,coords);
188 return computeVolSurfOfCell<ConnType,numPolConn,2>(type,connec,lgth,coords);
190 return computeVolSurfOfCell<ConnType,numPolConn,1>(type,connec,lgth,coords);
191 throw INTERP_KERNEL::Exception("Invalid spaceDim specified : must be 1, 2 or 3");
195 template<class ConnType, NumberingPolicy numPol,int SPACEDIM>
196 void computeBarycenter(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, double *res)
204 std::copy(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]),
205 coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]+1),res);
206 std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]),res,std::plus<double>());
207 std::transform(res,res+SPACEDIM,res,std::bind2nd(std::multiplies<double>(),0.5));
214 std::copy(coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]),
215 coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[0]+1),res);
216 std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[1]),res,std::plus<double>());
217 std::transform(res,res+SPACEDIM,coords+SPACEDIM*OTT<ConnType,numPol>::coo2C(connec[2]),res,std::plus<double>());
218 std::transform(res,res+SPACEDIM,res,std::bind2nd(std::multiplies<double>(),1./3.));
225 computePolygonBarycenter2D<ConnType,numPol>(connec,lgth,coords,res);
227 computePolygonBarycenter3D<ConnType,numPol>(connec,lgth,coords,res);
229 throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
235 computePolygonBarycenter2D<ConnType,numPol>(connec,lgth/2,coords,res);
237 computePolygonBarycenter3D<ConnType,numPol>(connec,lgth/2,coords,res);
239 throw INTERP_KERNEL::Exception("Impossible spacedim linked to cell 2D Cell !");
244 res[0]=coords[3*OTT<ConnType,numPol>::coo2C(connec[0])];
245 res[1]=coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+1];
246 res[2]=coords[3*OTT<ConnType,numPol>::coo2C(connec[0])+2];
247 res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[1])];
248 res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[1])+1];
249 res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[1])+2];
250 res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[2])];
251 res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[2])+1];
252 res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[2])+2];
253 res[0]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[3])];
254 res[1]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[3])+1];
255 res[2]+=coords[3*OTT<ConnType,numPol>::coo2C(connec[3])+2];
256 res[0]*=0.25; res[1]*=0.25; res[2]*=0.25;
262 computePolygonBarycenter3D<ConnType,numPol>(connec,lgth-1,coords,tmp);
263 res[0]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[4])]+3.*tmp[0])/4.;
264 res[1]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[4])+1]+3.*tmp[1])/4.;
265 res[2]=(coords[3*OTT<ConnType,numPol>::coo2C(connec[4])+2]+3.*tmp[2])/4.;
271 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,
272 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,
273 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,
274 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,
275 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,
276 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 barycenterOfPolyhedron<ConnType,numPol>(conn,29,coords,res);
284 OTT<ConnType,numPol>::coo2C(connec[0]),OTT<ConnType,numPol>::coo2C(connec[1]),OTT<ConnType,numPol>::coo2C(connec[2]),-1,
285 OTT<ConnType,numPol>::coo2C(connec[3]),OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[4]),-1,
286 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,
287 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,
288 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 barycenterOfPolyhedron<ConnType,numPol>(conn,22,coords,res);
293 case INTERP_KERNEL::NORM_HEXGP12:
295 const int connecHexa12[43]={
296 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,
297 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,
298 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,
299 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,
300 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,
301 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,
302 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,
303 OTT<ConnType,numPol>::coo2C(connec[5]),OTT<ConnType,numPol>::coo2C(connec[11]),OTT<ConnType,numPol>::coo2C(connec[6]),OTT<ConnType,numPol>::coo2C(connec[0])};
304 barycenterOfPolyhedron<ConnType,numPol>(connecHexa12,43,coords,res);
309 barycenterOfPolyhedron<ConnType,numPol>(connec,lgth,coords,res);
313 throw INTERP_KERNEL::Exception("Not recognized cell type to get Barycenter on it !");
317 template<class ConnType, NumberingPolicy numPolConn>
318 void computeBarycenter2(NormalizedCellType type, const ConnType *connec, int lgth, const double *coords, int spaceDim, double *res)
321 return computeBarycenter<ConnType,numPolConn,3>(type,connec,lgth,coords,res);
323 return computeBarycenter<ConnType,numPolConn,2>(type,connec,lgth,coords,res);
325 return computeBarycenter<ConnType,numPolConn,1>(type,connec,lgth,coords,res);
326 throw INTERP_KERNEL::Exception("Invalid spaceDim specified for compute barycenter : must be 1, 2 or 3");