1 // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege
2 //
3 // See the LICENSE.txt file for license information. Please report all
4 // issues on https://gitlab.onelab.info/getdp/getdp/issues.
5 //
6 // Contributor(s):
7 // Christophe Trophime
8 //
9
10 #include "ProData.h"
11 #include "Message.h"
12
13 #define SQU(a) ((a)*(a))
14
15 #define NoEdge Message::Error("Missing Edge Entity in Element %d", Element->Num)
16
17 /* ------------------------------------------------------------------------ */
18 /* B F _ E d g e */
19 /* ------------------------------------------------------------------------ */
20
21 #define WrongNumEdge Message::Error("Wrong Edge number in 'BF_Edge'")
22
BF_Edge(struct Element * Element,int NumEdge,double u,double v,double w,double s[])23 void BF_Edge(struct Element * Element, int NumEdge,
24 double u, double v, double w, double s[])
25 {
26 switch (Element->Type) {
27 case LINE : case LINE_2 :
28 case LINE_3 : case LINE_4 :
29 switch(NumEdge) {
30 case 1 : s[0] = 0.5 ; s[1] = 0. ; s[2] = 0. ; break ;
31 default : WrongNumEdge ;
32 }
33 break ;
34
35 case TRIANGLE : case TRIANGLE_2 :
36 case TRIANGLE_3 : case TRIANGLE_4 :
37 switch(NumEdge) {
38 case 1 : s[0] = 1.-v ; s[1] = u ; s[2] = 0. ; break ;
39 case 2 : s[0] = v ; s[1] = 1.-u ; s[2] = 0. ; break ;
40 case 3 : s[0] = -v ; s[1] = u ; s[2] = 0. ; break ;
41 default : WrongNumEdge ;
42 }
43 break ;
44
45 case QUADRANGLE : case QUADRANGLE_2 : case QUADRANGLE_2_8N :
46 case QUADRANGLE_3 : case QUADRANGLE_4 :
47 switch(NumEdge) {
48 case 1 : s[0] = 0.25 * (1.-v) ; s[1] = 0. ; s[2] = 0. ; break ;
49 case 2 : s[0] = 0. ; s[1] = 0.25 * (1.-u) ; s[2] = 0. ; break ;
50 case 3 : s[0] = 0. ; s[1] = 0.25 * (1.+u) ; s[2] = 0. ; break ;
51 case 4 : s[0] = -0.25 * (1.+v) ; s[1] = 0. ; s[2] = 0. ; break ;
52 default : WrongNumEdge ;
53 }
54 break ;
55
56 case TETRAHEDRON : case TETRAHEDRON_2 :
57 case TETRAHEDRON_3 : case TETRAHEDRON_4 :
58 switch(NumEdge) {
59 case 1 : s[0] = 1.-v-w ; s[1] = u ; s[2] = u ; break ;
60 case 2 : s[0] = v ; s[1] = 1.-u-w ; s[2] = v ; break ;
61 case 3 : s[0] = w ; s[1] = w ; s[2] = 1.-u-v ; break ;
62 case 4 : s[0] = -v ; s[1] = u ; s[2] = 0. ; break ;
63 case 5 : s[0] = -w ; s[1] = 0. ; s[2] = u ; break ;
64 case 6 : s[0] = 0. ; s[1] = -w ; s[2] = v ; break ;
65 default : WrongNumEdge ;
66 }
67 break ;
68
69 case HEXAHEDRON : case HEXAHEDRON_2 : case HEXAHEDRON_2_20N :
70 case HEXAHEDRON_3 : case HEXAHEDRON_4 :
71 switch(NumEdge) {
72 case 1 : s[0] = 0.125 * (1.-v) * (1.-w) ; s[1] = 0. ; s[2] = 0. ; break ;
73 case 6 : s[0] = -0.125 * (1.+v) * (1.-w) ; s[1] = 0. ; s[2] = 0. ; break ;
74 case 9 : s[0] = 0.125 * (1.-v) * (1.+w) ; s[1] = 0. ; s[2] = 0. ; break ;
75 case 12 : s[0] = -0.125 * (1.+v) * (1.+w) ; s[1] = 0. ; s[2] = 0. ; break ;
76
77 case 2 : s[0] = 0. ; s[1] = 0.125 * (1.-u) * (1.-w) ; s[2] = 0. ; break ;
78 case 4 : s[0] = 0. ; s[1] = 0.125 * (1.+u) * (1.-w) ; s[2] = 0. ; break ;
79 case 10 : s[0] = 0. ; s[1] = 0.125 * (1.-u) * (1.+w) ; s[2] = 0. ; break ;
80 case 11 : s[0] = 0. ; s[1] = 0.125 * (1.+u) * (1.+w) ; s[2] = 0. ; break ;
81
82 case 3 : s[0] = 0. ; s[1] = 0. ; s[2] = 0.125 * (1.-u) * (1.-v) ; break ;
83 case 5 : s[0] = 0. ; s[1] = 0. ; s[2] = 0.125 * (1.+u) * (1.-v) ; break ;
84 case 7 : s[0] = 0. ; s[1] = 0. ; s[2] = 0.125 * (1.+u) * (1.+v) ; break ;
85 case 8 : s[0] = 0. ; s[1] = 0. ; s[2] = 0.125 * (1.-u) * (1.+v) ; break ;
86 default : WrongNumEdge ;
87 }
88 break ;
89
90 case PRISM : case PRISM_2 : case PRISM_2_15N :
91 case PRISM_3 : case PRISM_4 :
92 switch(NumEdge) {
93 case 1 : s[0] = 0.5 * (1.-v) * (1.-w) ;
94 s[1] = 0.5 * u * (1.-w) ;
95 s[2] = 0. ; break ;
96 case 2 : s[0] = 0.5 * v * (1.-w) ;
97 s[1] = 0.5 * (1.-u) * (1.-w) ;
98 s[2] = 0. ; break ;
99 case 3 : s[0] = 0. ;
100 s[1] = 0. ;
101 s[2] = 0.5 * (1.-u-v) ; break ;
102 case 4 : s[0] = -0.5 * v * (1.-w) ;
103 s[1] = 0.5 * u * (1.-w) ;
104 s[2] = 0. ; break ;
105 case 5 : s[0] = 0. ;
106 s[1] = 0. ;
107 s[2] = 0.5 * u ; break ;
108 case 6 : s[0] = 0. ;
109 s[1] = 0. ;
110 s[2] = 0.5 * v ; break ;
111 case 7 : s[0] = 0.5 * (1.-v) * (1.+w) ;
112 s[1] = 0.5 * u * (1.+w) ;
113 s[2] = 0. ; break ;
114 case 8 : s[0] = 0.5 * v * (1.+w) ;
115 s[1] = 0.5 * (1.-u) * (1.+w) ;
116 s[2] = 0. ; break ;
117 case 9 : s[0] = -0.5 * v * (1.+w) ;
118 s[1] = 0.5 * u * (1.+w) ;
119 s[2] = 0. ; break ;
120 default : WrongNumEdge ;
121 }
122 break ;
123
124 case PYRAMID : case PYRAMID_2 : case PYRAMID_2_13N :
125 case PYRAMID_3 : // case PYRAMID_4
126 if (w != 1){
127 switch(NumEdge) {
128 case 1 : s[0] = 0.25 * (1 - v - w) ;
129 s[1] = 0. ;
130 s[2] = 0.25 * (u - u * v / (1. - w)) ; break ;
131 case 2 : s[0] = 0. ;
132 s[1] = 0.25 * (1 - u - w) ;
133 s[2] = 0.25 * (v - u * v / (1. - w)) ; break ;
134 case 4 : s[0] = 0. ;
135 s[1] = 0.25 * (1 + u - w) ;
136 s[2] = 0.25 * (v + u * v / (1. - w)) ; break ;
137 case 6 : s[0] = -0.25 * (1 + v - w) ;
138 s[1] = 0. ;
139 s[2] = -0.25 * (u + u * v / (1. - w)) ; break ;
140 case 3 : s[0] = 0.25 * (w - v * w / (1. - w)) ;
141 s[1] = 0.25 * (w - u * w / (1. - w)) ;
142 s[2] = 0.25 * (1. - u - v + u * v / SQU(1. - w) - 2 * u * v * w / SQU(1. - w)) ; break ;
143 case 5 : s[0] = -0.25 * (w - v * w / (1. - w)) ;
144 s[1] = 0.25 * (w + u * w / (1. - w)) ;
145 s[2] = 0.25 * (1. + u - v - u * v / SQU(1. - w) + 2 * u * v * w / SQU(1. - w)) ; break ;
146 case 7 : s[0] = -0.25 * (w + v * w / (1. - w)) ;
147 s[1] = -0.25 * (w + u * w / (1. - w)) ;
148 s[2] = 0.25 * (1. + u + v + u * v / SQU(1. - w) - 2 * u * v * w / SQU(1. - w)) ; break ;
149 case 8 : s[0] = 0.25 * (w + v * w / (1. - w)) ;
150 s[1] = -0.25 * (w - u * w / (1. - w)) ;
151 s[2] = 0.25 * (1. - u + v - u * v / SQU(1. - w) + 2 * u * v * w / SQU(1. - w)) ; break ;
152 default : WrongNumEdge ;
153 }
154 }
155 else
156 switch(NumEdge) {
157 case 1 : s[0] = -0.25 * v ;
158 s[1] = 0. ;
159 s[2] = 0.25 * u ; break ;
160 case 2 : s[0] = 0. ;
161 s[1] = -0.25 * u ;
162 s[2] = 0.25 * v ; break ;
163 case 4 : s[0] = 0. ;
164 s[1] = 0.25 * u ;
165 s[2] = 0.25 * v ; break ;
166 case 6 : s[0] = -0.25 * v ;
167 s[1] = 0. ;
168 s[2] = -0.25 * u ; break ;
169 case 3 : s[0] = 0.25 ;
170 s[1] = 0.25 ;
171 s[2] = 0.25 * (1. - u - v) ; break ;
172 case 5 : s[0] = -0.25 ;
173 s[1] = 0.25 ;
174 s[2] = 0.25 * (1. + u - v) ; break ;
175 case 7 : s[0] = -0.25 ;
176 s[1] = -0.25 ;
177 s[2] = 0.25 * (1. + u + v) ; break ;
178 case 8 : s[0] = 0.25 ;
179 s[1] = -0.25 ;
180 s[2] = 0.25 * (1. - u + v) ; break ;
181 default : WrongNumEdge ;
182 }
183 break ;
184
185 default :
186 Message::Error("Unknown type of Element in BF_Edge");
187 break;
188 }
189
190 if (!Element->GeoElement->NumEdges) NoEdge ;
191
192 if (Element->GeoElement->NumEdges[NumEdge-1] < 0) {
193 s[0] = - s[0] ; s[1] = - s[1] ; s[2] = - s[2] ;
194 }
195 }
196
197 #undef WrongNumEdge
198
199 /* ------------------------------------------------------------------------ */
200 /* B F _ C u r l E d g e */
201 /* ------------------------------------------------------------------------ */
202
203 #define WrongNumEdge Message::Error("Wrong Edge number in 'BF_CurlEdge'")
204
BF_CurlEdge(struct Element * Element,int NumEdge,double u,double v,double w,double s[])205 void BF_CurlEdge(struct Element * Element, int NumEdge,
206 double u, double v, double w, double s[])
207 {
208 switch (Element->Type) {
209 case LINE : case LINE_2 :
210 case LINE_3 : case LINE_4 :
211 switch(NumEdge) {
212 case 1 : s[0] = 0. ; s[1] = 0. ; s[2] = 0. ; break ;
213 default : WrongNumEdge ;
214 }
215 break ;
216
217 case TRIANGLE : case TRIANGLE_2 :
218 case TRIANGLE_3 : case TRIANGLE_4 :
219 switch(NumEdge) {
220 case 1 : s[0] = 0. ; s[1] = 0. ; s[2] = 2. ; break ;
221 case 2 : s[0] = 0. ; s[1] = 0. ; s[2] = -2. ; break ;
222 case 3 : s[0] = 0. ; s[1] = 0. ; s[2] = 2. ; break ;
223 default : WrongNumEdge ;
224 }
225 break ;
226
227 case QUADRANGLE : case QUADRANGLE_2 : case QUADRANGLE_2_8N :
228 case QUADRANGLE_3 : case QUADRANGLE_4 :
229 switch(NumEdge) {
230 case 1 : s[0] = 0. ; s[1] = 0. ; s[2] = 0.25 ; break ;
231 case 2 : s[0] = 0. ; s[1] = 0. ; s[2] = -0.25 ; break ;
232 case 3 : s[0] = 0. ; s[1] = 0. ; s[2] = 0.25 ; break ;
233 case 4 : s[0] = 0. ; s[1] = 0. ; s[2] = 0.25 ; break ;
234 default : WrongNumEdge ;
235 }
236 break ;
237
238 case TETRAHEDRON : case TETRAHEDRON_2 :
239 case TETRAHEDRON_3 : case TETRAHEDRON_4 :
240 switch(NumEdge) {
241 case 1 : s[0] = 0. ; s[1] = -2. ; s[2] = 2. ; break ;
242 case 2 : s[0] = 2. ; s[1] = 0. ; s[2] = -2. ; break ;
243 case 3 : s[0] = -2. ; s[1] = 2. ; s[2] = 0. ; break ;
244 case 4 : s[0] = 0. ; s[1] = 0. ; s[2] = 2. ; break ;
245 case 5 : s[0] = 0. ; s[1] = -2. ; s[2] = 0. ; break ;
246 case 6 : s[0] = 2. ; s[1] = 0. ; s[2] = 0. ; break ;
247 default : WrongNumEdge ;
248 }
249 break ;
250
251 case HEXAHEDRON : case HEXAHEDRON_2 : case HEXAHEDRON_2_20N :
252 case HEXAHEDRON_3 : case HEXAHEDRON_4 :
253 switch(NumEdge) {
254 case 1 : s[0] = 0. ; s[1] = 0.125*(v-1.) ; s[2] = 0.125*(1.-w) ; break ;
255 case 6 : s[0] = 0. ; s[1] = 0.125*(v+1.) ; s[2] = 0.125*(1.-w) ; break ;
256 case 9 : s[0] = 0. ; s[1] = 0.125*(1.-v) ; s[2] = 0.125*(1.+w) ; break ;
257 case 12 : s[0] = 0. ; s[1] =-0.125*(v+1.) ; s[2] = 0.125*(1.+w) ; break ;
258
259 case 2 : s[0] = 0.125*(1.-u) ; s[1] = 0. ; s[2] = 0.125*(w-1.) ; break ;
260 case 4 : s[0] = 0.125*(1.+u) ; s[1] = 0. ; s[2] = 0.125*(1.-w) ; break ;
261 case 10 : s[0] = 0.125*(u-1.) ; s[1] = 0. ; s[2] =-0.125*(w+1.) ; break ;
262 case 11 : s[0] =-0.125*(1.+u) ; s[1] = 0. ; s[2] = 0.125*(w+1.) ; break ;
263
264 case 3 : s[0] = 0.125*(u-1.) ; s[1] = 0.125*(1.-v) ; s[2] = 0. ; break ;
265 case 5 : s[0] =-0.125*(u+1.) ; s[1] = 0.125*(v-1.) ; s[2] = 0. ; break ;
266 case 7 : s[0] = 0.125*(u+1.) ; s[1] =-0.125*(1.+v) ; s[2] = 0. ; break ;
267 case 8 : s[0] = 0.125*(1.-u) ; s[1] = 0.125*(1.+v) ; s[2] = 0. ; break ;
268 default : WrongNumEdge ;
269 }
270 break ;
271
272 case PRISM : case PRISM_2 : case PRISM_2_15N :
273 case PRISM_3 : case PRISM_4 :
274 switch(NumEdge) {
275 case 1 : s[0] = 0.5*u ; s[1] = 0.5*(v-1.) ; s[2] = 1.-w ; break ;
276 case 2 : s[0] = 0.5*(1.-u) ; s[1] = -0.5*v ; s[2] = w-1. ; break ;
277 case 3 : s[0] = -0.5 ; s[1] = 0.5 ; s[2] = 0. ; break ;
278
279 case 4 : s[0] = 0.5*u ; s[1] = 0.5*v ; s[2] = 1.-w ; break ;
280 case 5 : s[0] = 0. ; s[1] = -0.5 ; s[2] = 0. ; break ;
281 case 6 : s[0] = 0.5 ; s[1] = 0. ; s[2] = 0. ; break ;
282
283 case 7 : s[0] = -0.5*u ; s[1] = 0.5*(1.-v) ; s[2] = 1.+w ; break ;
284 case 8 : s[0] = 0.5*(u-1.) ; s[1] = 0.5*v ; s[2] = -1.-w ; break ;
285 case 9 : s[0] = -0.5*u ; s[1] = -0.5*v ; s[2] = 1.+w ; break ;
286 default : WrongNumEdge ;
287 }
288 break ;
289
290 case PYRAMID : case PYRAMID_2 : case PYRAMID_2_13N :
291 case PYRAMID_3 : // case PYRAMID_4
292 if (w != 1){
293 switch(NumEdge) {
294 case 1 : s[0] = -0.25 * u / (1. - w) ; s[1] = -0.5 + 0.25 * v / (1. - w) ; s[2] = 0.25 ; break ;
295 case 2 : s[0] = 0.5 - 0.25 * u / (1. - w) ; s[1] = 0.25 * v / (1. - w) ; s[2] = -0.25 ; break ;
296 case 4 : s[0] = 0.5 + 0.25 * u / (1. - w) ; s[1] = -0.25 * v / (1. - w) ; s[2] = 0.25 ; break ;
297 case 6 : s[0] = -0.25 * u / (1. - w) ; s[1] = 0.5 + 0.25 * v / (1. - w) ; s[2] = 0.25 ; break ;
298 case 3 : s[0] = -0.5 * (1. - u / (1. - w)) ; s[1] = 0.5 * (1. - v / (1. - w)) ; s[2] = 0. ; break;
299 case 5 : s[0] = -0.5 * (1. + u / (1. - w)) ; s[1] = -0.5 * (1. - v / (1. - w)) ; s[2] = 0. ; break;
300 case 7 : s[0] = 0.5 * (1. + u / (1. - w)) ; s[1] = -0.5 * (1. + v / (1. - w)) ; s[2] = 0. ; break;
301 case 8 : s[0] = 0.5 * (1. - u / (1. - w)) ; s[1] = 0.5 * (1. + v / (1. - w)) ; s[2] = 0. ; break;
302 default : WrongNumEdge ;
303 }
304 }
305 else {
306 switch(NumEdge) {
307 case 1 : s[0] = 0. ; s[1] = -0.5 ; s[2] = 0.25 ; break ;
308 case 2 : s[0] = 0.5 ; s[1] = 0. ; s[2] = -0.25 ; break ;
309 case 4 : s[0] = 0.5 ; s[1] = 0. ; s[2] = 0.25 ; break ;
310 case 6 : s[0] = 0. ; s[1] = 0.5 ; s[2] = 0.25 ; break ;
311 case 3 : s[0] = -0.5 ; s[1] = 0.5 ; s[2] = 0. ; break;
312 case 5 : s[0] = -0.5 ; s[1] = -0.5 ; s[2] = 0. ; break;
313 case 7 : s[0] = 0.5 ; s[1] = -0.5 ; s[2] = 0. ; break;
314 case 8 : s[0] = 0.5 ; s[1] = 0.5 ; s[2] = 0. ; break;
315 default : WrongNumEdge ;
316 }
317 }
318 break ;
319
320
321 default :
322 Message::Error("Unknown type of Element in BF_CurlEdge");
323 break;
324 }
325
326 if (!Element->GeoElement->NumEdges) NoEdge ;
327
328 if (Element->GeoElement->NumEdges[NumEdge-1] < 0) {
329 s[0] = - s[0] ; s[1] = - s[1] ; s[2] = - s[2] ;
330 }
331 }
332
333 #undef WrongNumEdge
334
335 #undef NoEdge
336