1 /*
2 XLiFE++ is an extended library of finite elements written in C++
3 Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
4
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13 You should have received a copy of the GNU General Public License
14 along with this program. If not, see <http://www.gnu.org/licenses/>.
15 */
16
17 /*!
18 \file GeomRefPrism.cpp
19 \authors D. Martin, N.Kielbasiewicz
20 \since 15 dec 2002
21 \date 19 sep 2012
22
23 \brief Implementation of xlifepp::GeomRefPrism class members and related functions
24 */
25
26 #include "GeomRefPrism.hpp"
27 #include "utils.h"
28
29 namespace xlifepp
30 {
31
32 //--------------------------------------------------------------------------------
33 // GeomRefPrism constructor for Geometric Reference Element
34 //--------------------------------------------------------------------------------
GeomRefPrism()35 GeomRefPrism::GeomRefPrism()
36 : GeomRefElement(_prism, 0.5, over3_, 6, 9)
37 //! using 3d base constructor with shape, volume, centroid coords, number of vertices, edges
38 {
39 trace_p->push("GeomRefPrism::GeomRefPrism");
40 centroid_[2] = 0.5;
41 // coordinates of vertices
42 std::vector<real_t>::iterator it_v(vertices_.begin());
43 vertex(it_v, 1., 0., 0.);
44 vertex(it_v, 0., 1., 0.);
45 vertex(it_v, 0., 0., 0.);
46 vertex(it_v, 1., 0., 1.);
47 vertex(it_v, 0., 1., 1.);
48 vertex(it_v, 0., 0., 1.);
49 // vertex numbering on edges
50 sideOfSideNumbering();
51 // vertex numbering and oriented edge numbering on faces
52 sideNumbering();
53
54 trace_p->pop();
55 }
56
57 /*
58 --------------------------------------------------------------------------------
59 sideOfSideNumbering
60 defines Geometric Reference Element local numbering of edge vertices
61 sideNumbering
62 defines Geometric Reference Element local numbering of vertices on faces
63 and of edges on faces
64 --------------------------------------------------------------------------------
65 */
sideOfSideNumbering()66 void GeomRefPrism::sideOfSideNumbering()
67 {
68 sideOfSideVertexNumbers_[0].push_back(1);
69 sideOfSideVertexNumbers_[0].push_back(2);
70 sideOfSideVertexNumbers_[1].push_back(2);
71 sideOfSideVertexNumbers_[1].push_back(3);
72 sideOfSideVertexNumbers_[2].push_back(3);
73 sideOfSideVertexNumbers_[2].push_back(1);
74 sideOfSideVertexNumbers_[3].push_back(4);
75 sideOfSideVertexNumbers_[3].push_back(5);
76 sideOfSideVertexNumbers_[4].push_back(5);
77 sideOfSideVertexNumbers_[4].push_back(6);
78 sideOfSideVertexNumbers_[5].push_back(6);
79 sideOfSideVertexNumbers_[5].push_back(4);
80 sideOfSideVertexNumbers_[6].push_back(1);
81 sideOfSideVertexNumbers_[6].push_back(4);
82 sideOfSideVertexNumbers_[7].push_back(2);
83 sideOfSideVertexNumbers_[7].push_back(5);
84 sideOfSideVertexNumbers_[8].push_back(3);
85 sideOfSideVertexNumbers_[8].push_back(6);
86
87 sideOfSideNumbers_[0].push_back(1);
88 sideOfSideNumbers_[0].push_back(-3);
89 sideOfSideNumbers_[0].push_back(-2);
90 sideOfSideNumbers_[1].push_back(1);
91 sideOfSideNumbers_[1].push_back(7);
92 sideOfSideNumbers_[1].push_back(4);
93 sideOfSideNumbers_[1].push_back(-8);
94 sideOfSideNumbers_[2].push_back(-9);
95 sideOfSideNumbers_[2].push_back(-5);
96 sideOfSideNumbers_[2].push_back(-8);
97 sideOfSideNumbers_[2].push_back(2);
98 sideOfSideNumbers_[3].push_back(-3);
99 sideOfSideNumbers_[3].push_back(-7);
100 sideOfSideNumbers_[3].push_back(6);
101 sideOfSideNumbers_[3].push_back(9);
102 sideOfSideNumbers_[4].push_back(4);
103 sideOfSideNumbers_[4].push_back(-5);
104 sideOfSideNumbers_[4].push_back(-6);
105
106 }
107
sideNumbering()108 void GeomRefPrism::sideNumbering()
109 {
110 for (number_t i = 1; i < nbSides_ - 1; i++)
111 {
112 sideShapeTypes_[i] = _quadrangle;
113 }
114 sideShapeTypes_[0] = _triangle;
115 sideShapeTypes_[nbSides_ - 1] = _triangle;
116
117 sideVertexNumbers_[0].push_back(1);
118 sideVertexNumbers_[0].push_back(3);
119 sideVertexNumbers_[0].push_back(2);
120 sideVertexNumbers_[1].push_back(1);
121 sideVertexNumbers_[1].push_back(2);
122 sideVertexNumbers_[1].push_back(5);
123 sideVertexNumbers_[1].push_back(4);
124 sideVertexNumbers_[2].push_back(2);
125 sideVertexNumbers_[2].push_back(3);
126 sideVertexNumbers_[2].push_back(6);
127 sideVertexNumbers_[2].push_back(5);
128 sideVertexNumbers_[3].push_back(3);
129 sideVertexNumbers_[3].push_back(1);
130 sideVertexNumbers_[3].push_back(4);
131 sideVertexNumbers_[3].push_back(6);
132 sideVertexNumbers_[4].push_back(4);
133 sideVertexNumbers_[4].push_back(5);
134 sideVertexNumbers_[4].push_back(6);
135
136 }
137
138 //! Returns edge length, face area or element volume
measure(const dimen_t d,const number_t sideNum) const139 real_t GeomRefPrism::measure(const dimen_t d, const number_t sideNum) const
140 {
141 real_t ms = 0.;
142 switch (d)
143 {
144 case 0: ms = 1.; break;
145 case 1:
146 switch (sideNum)
147 {
148 case 1: case 4: ms = sqrtOf2_; break;
149 case 2: case 3: case 5: case 6: case 7: case 8: case 9: ms = 1.; break;
150 default: noSuchSideOfSide(sideNum);
151 }
152 break;
153 case 2:
154 switch (sideNum)
155 {
156 case 1: case 5: ms = .5; break;
157 case 3: case 4: ms = 1.; break;
158 case 2: ms = sqrtOf2_; break;
159 default: noSuchSide(sideNum);
160 }
161 break;
162 case 3: ms = measure_; break;
163 default: break;
164 }
165 return ms;
166 }
167
168 // /*
169 // --------------------------------------------------------------------------------
170 // Returns a quadrature rule built on an edge from a 1d quadrature formula
171 // or a quadrature rule built on an face from a triangle or quadrangle quadrature formula
172 // --------------------------------------------------------------------------------
173 // */
174 // void GeomRefPrism::sideQuadrature(const QuadratureRule& q1, QuadratureRule& qr, const number_t sideNum, const dimen_t d) const
175 // {
176 // /*
177 // q1 : input 1d or 2d (triangle or quadrangle) quadrature rule
178 // q2 : ouput 3d quadrature rule set on the edge or face
179 // sideNum : local number of edge (sideNum = 1,2,...,9) or face (sideNum = 1,2,3,4,5)
180 // d : side dimension (d = 1 for an edge, d = 2 for a face)
181 // */
182 // std::vector<real_t>::const_iterator c_1(q1.point()), w_1(q1.weight());
183 // std::vector<real_t>::iterator c_i(qr.point()), w_i(qr.weight());
184
185 // switch ( d )
186 // {
187 // case 2:
188 // switch ( sideNum )
189 // {
190 // case 1: // face (x3 = 0)
191 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1 += 2 )
192 // qr.point(c_i, *c_1, *(c_1+1), 0., w_i, *w_1);
193 // break;
194 // case 2: // face (x1 + x2 = 1)
195 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1 += 2 )
196 // qr.point(c_i, *c_1, 1.-*c_1, *(c_1+1), w_i, *w_1*sqrt_of_2_);
197 // break;
198 // case 3: // face (x1 = 0)
199 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1 += 2 )
200 // qr.point(c_i, 0., *c_1, *(c_1+1), w_i, *w_1);
201 // break;
202 // case 4: // face (x2 = 0)
203 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1 += 2 )
204 // qr.point(c_i, *c_1, 0., *(c_1+1), w_i, *w_1);
205 // break;
206 // case 5: // face (x3 = 1)
207 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1 += 2 )
208 // qr.point(c_i, *c_1, *(c_1+1), 1., w_i, *w_1);
209 // break;
210 // default: noSuchFace(sideNum); break;
211 // }
212 // break;
213 // case 1:
214 // switch ( sideNum )
215 // {
216 // // edges in plane x3 = 0
217 // case 1: // edge (x1 + x2 = 1 & x3 = 0)
218 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
219 // qr.point(c_i, *c_1, 1.-*c_1, 0., w_i, *w_1*sqrt_of_2_);
220 // break;
221 // case 2: // edge (x1 = 0 & x3 = 0)
222 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
223 // qr.point(c_i, 0., *c_1, 0., w_i, *w_1);
224 // break;
225 // case 3: // edge (x2 = 0 & x3 = 0)
226 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
227 // qr.point(c_i, *c_1, 0., 0., w_i, *w_1);
228 // break;
229 // // edges in plane x3 = 1
230 // case 4: // edge (x1 + x2 = 1 & x3 = 1)
231 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
232 // qr.point(c_i, *c_1, 1.-*c_1, 1., w_i, *w_1*sqrt_of_2_);
233 // break;
234 // case 5: // edge (x1 = 0 & x3 = 1)
235 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
236 // qr.point(c_i, 0., *c_1, 1., w_i, *w_1);
237 // break;
238 // case 6: // edge (x2 = 0 & x3 = 1)
239 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
240 // qr.point(c_i, *c_1, 0., 1., w_i, *w_1);
241 // break;
242 // // edges orthogonal to x3 = 0 plane
243 // case 7: // edge (x1 = 1 & x2 = 0)
244 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
245 // qr.point(c_i, 1., 0., *c_1, w_i, *w_1);
246 // break;
247 // case 8: // edge (x1 = 0 & x2 = 1)
248 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
249 // qr.point(c_i, 0., 1., *c_1, w_i, *w_1);
250 // break;
251 // case 9: // edge (x1 = 0 & x2 = 0)
252 // for ( size_t i = 0; i < q1.size(); i++, w_1++, c_1++ )
253 // qr.point(c_i, 0., 0., *c_1, w_i, *w_1);
254 // break;
255 // default: noSuchEdge(sideNum);
256 // break;
257 // }
258 // default: break;
259 // }
260 // }
261 //! returns tangent vector(s) on edge sideNum (=1,2,...,9) or face sideNum (=1,2,3,4,5)
tangentVector(const std::vector<real_t> & jacobianMatrix,std::vector<std::vector<real_t>> & tgv,const number_t sideNum,const dimen_t d) const262 void GeomRefPrism::tangentVector(const std::vector<real_t>& jacobianMatrix, std::vector <std::vector<real_t> >& tgv, const number_t sideNum, const dimen_t d) const
263 {
264 std::vector<real_t>::iterator it_tgv0(tgv[0].begin()), it_tgv1(it_tgv0);
265 if ( d > 1) { it_tgv1 = tgv[1].begin(); }
266 // std::vector<real_t>::const_iterator it_jm(jacobianMatrix.begin());
267
268 noSuchFunction("tangentVector");
269
270 // switch ( d )
271 // {
272 // case 2:
273 // switch ( sideNum )
274 // {
275 // case 1: //
276 // break;
277 // case 2: //
278 // break;
279 // case 3: //
280 // break;
281 // case 4: //
282 // break;
283 // case 5: //
284 // break;
285 // default: noSuchSide(sideNum); break;
286 // }
287 // case 1:
288 // switch ( sideNum )
289 // {
290 // case 1: //
291 // break;
292 // case 2: //
293 // break;
294 // case 3: //
295 // break;
296 // case 4: //
297 // break;
298 // case 5: //
299 // break;
300 // case 6: //
301 // break;
302 // case 7: //
303 // break;
304 // case 8: //
305 // break;
306 // case 9: //
307 // break;
308 // default: noSuchSideOfSide(sideNum); break;
309 // }
310 // break;
311 // default:
312 // break;
313 // }
314 }
315
316 //! returns local number of edge bearing vertices with local numbers v1 and v2
sideWithVertices(const number_t vn1,const number_t vn2) const317 number_t GeomRefPrism::sideWithVertices(const number_t vn1, const number_t vn2) const
318 {
319 if(vn1==vn2) noSuchSide(vn1,vn2);
320 number_t v1=vn1, v2=vn2;
321 if(v1>v2) {v1=vn2; v2=vn1;}
322 switch(v1)
323 {
324 case 1 :
325 {
326 if(v2==2) return 1;
327 if(v2==3) return 3;
328 if(v2==4) return 7;
329 }
330 break;
331 case 2 :
332 {
333 if(v2==3) return 2;
334 if(v2==5) return 8;
335 }
336 break;
337 case 3 : if(v2==6) return 9; break;
338 case 4 :
339 {
340 if(v2==5) return 4;
341 if(v2==6) return 6;
342 }
343 break;
344 case 5 : if(v2==6) return 5;
345 }
346 noSuchSide(vn1,vn2);
347 return 0;
348 }
349 /*! returns local number of face bearing vertices with local numbers v1, v2, v3 and v4 if not null
350 1: 1 2 3 100v1+10v2+v3 123
351 2: 1 2 5 4 124 125 145 245
352 3: 2 3 6 5 235 236 256 356
353 4: 1 3 6 4 134 136 146 346
354 5: 4 5 6 456
355 */
sideWithVertices(const number_t vn1,const number_t vn2,const number_t vn3,const number_t vn4) const356 number_t GeomRefPrism::sideWithVertices(const number_t vn1, const number_t vn2, const number_t vn3, const number_t vn4) const
357 {
358 std::set<number_t> vs; vs.insert(vn1); vs.insert(vn2); vs.insert(vn3);
359 if(vn4>0)
360 {
361 vs.insert(vn4);
362 if(vs.size()!=4 || *vs.begin()<1 || *vs.rbegin()>6) noSuchSide(vn1,vn2,vn3,vn4);
363 }
364 else if(vs.size()!=3 || *vs.begin()<1 || *vs.rbegin()>6) noSuchSide(vn1,vn2,vn3);
365 std::map<number_t,number_t> hm;
366 std::map<number_t,number_t>::iterator itm;
367 number_t s=0;
368 std::set<number_t>::iterator itv=vs.begin();
369 if(vn4==0)
370 {
371 hm[123]=1; hm[124]=2; hm[125]=2; hm[145]=2; hm[245]=2; hm[235]=3; hm[236]=3; hm[256]=3; hm[356]=3;
372 hm[134]=4; hm[136]=4; hm[146]=4; hm[346]=4; hm[456]=5;
373 s=100**itv++; s+=10**itv++; s+=*itv;
374 itm=hm.find(s);
375 if(itm!=hm.end()) return itm->second; else noSuchSide(vn1,vn2,vn3);
376 }
377 hm[1245]=2; hm[2356]=3; hm[1346]=4;
378 s=1000**itv++; s+=100**itv++; s+=10**itv++;s+=*itv;
379 itm=hm.find(s);
380 if(itm!=hm.end()) return itm->second; else noSuchSide(vn1,vn2,vn3,vn4);
381 return 0;
382 }
383
384 //! test if a point belongs to current element
contains(std::vector<real_t> & p,real_t tol) const385 bool GeomRefPrism::contains(std::vector<real_t>& p, real_t tol) const
386 {
387 real_t x=p[0], y=p[1], z=p[2];
388 return (x >= -tol) && (x <= 1.+tol) && (y >= -tol) && (y <= 1.+tol) && (z >= -tol) && (z <= 1.+tol)
389 && (x+y<=1.+tol);
390 }
391
392 } // end of namespace xlifepp
393