BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
cl_IF_TET20.hpp
Go to the documentation of this file.
1/*
2 * BELFEM -- The Berkeley Lab Finite Element Framework
3 * Copyright (c) 2026, The Regents of the University of California,
4 * through Lawrence Berkeley National Laboratory (subject to receipt of any required
5 * approvals from the U.S. Dept. of Energy). All rights reserved.
6 *
7 * Developers: Christian Messe, Gregory Giard
8 *
9 * See the top-level LICENSE file for the complete license and disclaimer.
10 */
11
12#ifndef BELFEM_CL_IF_TET20_HPP
13#define BELFEM_CL_IF_TET20_HPP
14
15namespace belfem
16{
17 namespace fem
18 {
19//------------------------------------------------------------------------------
20
21 template<>
29
30//------------------------------------------------------------------------------
31
32 template<>
37 {
38 return ElementType::TET20;
39 }
40
41//------------------------------------------------------------------------------
42
43 template<>
44 void
48 {
49 aXiHat.set_size( 3, 20 );
50
51 const real a = 1.0/3.0;
52 const real b = 2.0/3.0;
53
54 aXiHat( 0, 0 ) = 1.0;
55 aXiHat( 1, 0 ) = 0.0;
56 aXiHat( 2, 0 ) = 0.0;
57
58 aXiHat( 0, 1 ) = 0.0;
59 aXiHat( 1, 1 ) = 0.0;
60 aXiHat( 2, 1 ) = 1.0;
61
62 aXiHat( 0, 2 ) = 0.0;
63 aXiHat( 1, 2 ) = 1.0;
64 aXiHat( 2, 2 ) = 0.0;
65
66 aXiHat( 0, 3 ) = 0.0;
67 aXiHat( 1, 3 ) = 0.0;
68 aXiHat( 2, 3 ) = 0.0;
69
70 aXiHat( 0, 4 ) = b;
71 aXiHat( 1, 4 ) = 0.0;
72 aXiHat( 2, 4 ) = a;
73
74 aXiHat( 0, 5 ) = a;
75 aXiHat( 1, 5 ) = 0.0;
76 aXiHat( 2, 5 ) = b;
77
78 aXiHat( 0, 6 ) = 0.0;
79 aXiHat( 1, 6 ) = a;
80 aXiHat( 2, 6 ) = b;
81
82 aXiHat( 0, 7 ) = 0.0;
83 aXiHat( 1, 7 ) = b;
84 aXiHat( 2, 7 ) = a;
85
86 aXiHat( 0, 8 ) = a;
87 aXiHat( 1, 8 ) = b;
88 aXiHat( 2, 8 ) = 0.0;
89
90 aXiHat( 0, 9 ) = b;
91 aXiHat( 1, 9 ) = a;
92 aXiHat( 2, 9 ) = 0.0;
93
94 aXiHat( 0, 10 ) = a;
95 aXiHat( 1, 10 ) = 0.0;
96 aXiHat( 2, 10 ) = 0.0;
97
98 aXiHat( 0, 11 ) = b;
99 aXiHat( 1, 11 ) = 0.0;
100 aXiHat( 2, 11 ) = 0.0;
101
102 aXiHat( 0, 12 ) = 0.0;
103 aXiHat( 1, 12 ) = a;
104 aXiHat( 2, 12 ) = 0.0;
105
106 aXiHat( 0, 13 ) = 0.0;
107 aXiHat( 1, 13 ) = b;
108 aXiHat( 2, 13 ) = 0.0;
109
110 aXiHat( 0, 14 ) = 0.0;
111 aXiHat( 1, 14 ) = 0.0;
112 aXiHat( 2, 14 ) = a;
113
114 aXiHat( 0, 15 ) = 0.0;
115 aXiHat( 1, 15 ) = 0.0;
116 aXiHat( 2, 15 ) = b;
117
118 aXiHat( 0, 16 ) = a;
119 aXiHat( 1, 16 ) = a;
120 aXiHat( 2, 16 ) = a;
121
122 aXiHat( 0, 17 ) = a;
123 aXiHat( 1, 17 ) = 0.0;
124 aXiHat( 2, 17 ) = a;
125
126 aXiHat( 0, 18 ) = a;
127 aXiHat( 1, 18 ) = a;
128 aXiHat( 2, 18 ) = 0.0;
129
130 aXiHat( 0, 19 ) = 0.0;
131 aXiHat( 1, 19 ) = a;
132 aXiHat( 2, 19 ) = a;
133 }
134//------------------------------------------------------------------------------
135
136 template<>
137 void
140 const Vector< real > & aXi,
141 Matrix< real > & aN ) const
142 {
143 const real xi = aXi( 0 );
144 const real eta = aXi( 1 );
145 const real zeta = aXi( 2 ) ;
146
147 const real tau = 1.0 - xi - eta - zeta;
148 const real psi = 3.0*(xi+eta+zeta);
149
150 aN.set_size( 1, 20 );
151
152 aN( 0, 0 ) = (xi*(9.0*xi*(xi-1.0)+ 2.0))* 0.5;
153 aN( 0, 1 ) = (zeta*(9.0*zeta*(zeta-1.0) + 2.0)) * 0.5;
154 aN( 0, 2 ) = (eta*(9.0*eta*(eta-1.0)+2.0)) * 0.5;
155 aN( 0, 3 ) = 0.5 * tau*( psi-2.0 )*( psi-1.0 );
156 aN( 0, 4 ) = 4.5*xi*zeta*(3.0*xi-1.0);
157 aN( 0, 5 ) = 4.5*xi*zeta*(3.0*zeta-1.0);
158 aN( 0, 6 ) = 4.5*eta*zeta*(3.0*zeta-1.0);
159 aN( 0, 7 ) = 4.5*eta*zeta*(3.0*eta-1.0);
160 aN( 0, 8 ) = 4.5*eta*xi*(3.0*eta-1.0);
161 aN( 0, 9 ) = 4.5*eta*xi*(3.0*xi-1.0);
162 aN( 0, 10 ) = 4.5*xi*tau*( 2.0-psi );
163 aN( 0, 11 ) = 4.5*xi*(3.0*xi-1.0)*tau;
164 aN( 0, 12 ) = 4.5*eta*tau*(2.0-psi);
165 aN( 0, 13 ) = 4.5*eta*(3.0*eta-1.0)*tau;
166 aN( 0, 14 ) = 4.5*zeta*tau*(2.0-psi);
167 aN( 0, 15 ) = 4.5*zeta*(3.0*zeta-1.0)*tau;
168 aN( 0, 16 ) = 27.0*xi*eta*zeta;
169 aN( 0, 17 ) = 27.0*xi*zeta*tau;
170 aN( 0, 18 ) = 27.0*xi*eta*tau;
171 aN( 0, 19 ) = 27.0*eta*zeta*tau;
172
173 }
174
175//------------------------------------------------------------------------------
176
177 template<>
178 void
181 const Vector< real > & aXi,
182 Matrix< real > & adNdXi ) const
183 {
184 const real xi = aXi( 0 );
185 const real eta = aXi( 1 );
186 const real zeta = aXi( 2 );
187
188 const real xi2 = xi*xi;
189 const real eta2 = eta*eta;
190
191 adNdXi.set_size( 3, 20 );
192
193 adNdXi( 0, 0 ) = 1.0+xi*(13.5*xi-9.0);
194 adNdXi( 1, 0 ) = 0.0;
195 adNdXi( 2, 0 ) = 0.0;
196
197 adNdXi( 0, 1 ) = 0.0;
198 adNdXi( 1, 1 ) = 0.0;
199 adNdXi( 2, 1 ) = 1.0+zeta*(13.5*zeta-9.0);
200
201 adNdXi( 0, 2 ) = 0.0;
202 adNdXi( 1, 2 ) = 1.0+eta*(13.5*eta-9.0);
203 adNdXi( 2, 2 ) = 0.0;
204
205 adNdXi( 0, 3 ) = xi*(18.0-27.0*zeta)-13.5*(xi2+eta2)-5.5+eta*(18.0-27.0*(xi+zeta))+zeta*(18.0-13.5*zeta);
206 adNdXi( 1, 3 ) = xi*(18.0-27.0*zeta)-13.5*(xi2+eta2)-5.5+eta*(18.0-27.0*(xi+zeta))+zeta*(18.0-13.5*zeta);
207 adNdXi( 2, 3 ) = xi*(18.0-27.0*zeta)-13.5*(xi2+eta2)-5.5+eta*(18.0-27.0*(xi+zeta))+zeta*(18.0-13.5*zeta);
208
209 adNdXi( 0, 4 ) = zeta*(27.0*xi-4.5);
210 adNdXi( 1, 4 ) = 0.0;
211 adNdXi( 2, 4 ) = xi*(13.5*xi-4.5);
212
213 adNdXi( 0, 5 ) = zeta*(13.5*zeta-4.5);
214 adNdXi( 1, 5 ) = 0.0;
215 adNdXi( 2, 5 ) = xi*(27.0*zeta-4.5);
216
217 adNdXi( 0, 6 ) = 0.0;
218 adNdXi( 1, 6 ) = zeta*(13.5*zeta-4.5);
219 adNdXi( 2, 6 ) = eta*(27.0*zeta-4.5);
220
221 adNdXi( 0, 7 ) = 0.0;
222 adNdXi( 1, 7 ) = zeta*(27.0*eta-4.5);
223 adNdXi( 2, 7 ) = eta*(13.5*eta-4.5);
224
225 adNdXi( 0, 8 ) = eta*(13.5*eta-4.5);
226 adNdXi( 1, 8 ) = xi*(27.0*eta-4.5);
227 adNdXi( 2, 8 ) = 0.0;
228
229 adNdXi( 0, 9 ) = eta*(27.0*xi-4.5);
230 adNdXi( 1, 9 ) = xi*(13.5*xi-4.5);
231 adNdXi( 2, 9 ) = 0.0;
232
233 adNdXi( 0, 10 ) = 9.0+13.5*eta2+40.5*xi2+zeta*(13.5*zeta-22.5)+eta*(54.0*xi+27.0*zeta-22.5)+xi*(54.0*zeta-45.0);
234 adNdXi( 1, 10 ) = xi*(27.0*(xi+eta+zeta)-22.5);
235 adNdXi( 2, 10 ) = xi*(27.0*(xi+eta+zeta)-22.5);
236
237 adNdXi( 0, 11 ) = xi*(36.0-40.5*xi-27.0*zeta)+eta*(4.5-27.0*xi)+4.5*zeta-4.5;
238 adNdXi( 1, 11 ) = (4.5-13.5*xi)*xi;
239 adNdXi( 2, 11 ) = (4.5-13.5*xi)*xi;
240
241 adNdXi( 0, 12 ) = eta*(27.0*(xi+eta+zeta)-22.5);
242 adNdXi( 1, 12 ) = 9.0+40.5*eta2+13.5*xi2+xi*(27.0*zeta-22.5)+eta*(54.0*xi+54.0*zeta-45.0)+zeta*(13.5*zeta-22.5);
243 adNdXi( 2, 12 ) = eta*(27.0*(xi+eta+zeta)-22.5);
244
245 adNdXi( 0, 13 ) = eta*(4.5-13.5*eta);
246 adNdXi( 1, 13 ) = 4.5*(xi+zeta)+eta*(36.0-40.5*eta-27.0*(xi+zeta))-4.5;
247 adNdXi( 2, 13 ) = eta*(4.5-13.5*eta);
248
249 adNdXi( 0, 14 ) = zeta*(27.0*(xi+eta+zeta)-22.5);
250 adNdXi( 1, 14 ) = zeta*(27.0*(xi+eta+zeta)-22.5);
251 adNdXi( 2, 14 ) = 9.0+13.5*(xi2+eta2)+zeta*(40.5*zeta-45.0)+xi*(54.0*zeta-22.5)+eta*(27.0*xi+54.0*zeta-22.5);
252
253 adNdXi( 0, 15 ) = zeta*(4.5-13.5*zeta);
254 adNdXi( 1, 15 ) = zeta*(4.5-13.5*zeta);
255 adNdXi( 2, 15 ) = xi*(4.5-27.0*zeta)+eta*(4.5-27.0*zeta)+zeta*(36.0-40.5*zeta)-4.5;
256
257 adNdXi( 0, 16 ) = 27.0*eta*zeta;
258 adNdXi( 1, 16 ) = 27.0*xi*zeta;
259 adNdXi( 2, 16 ) = 27.0*xi*eta;
260
261 adNdXi( 0, 17 ) = 27.0*zeta*(1.0-(2.0*xi+eta+zeta));
262 adNdXi( 1, 17 ) = -27.0*xi*zeta;
263 adNdXi( 2, 17 ) = 27.0*xi*(1.0-(xi+eta+2.0*zeta));
264
265 adNdXi( 0, 18 ) = 27.0*eta*(1.0 - (eta+zeta+2.0*xi));
266 adNdXi( 1, 18 ) = 27.0*xi*(1.0-(xi+2.0*eta+zeta));
267 adNdXi( 2, 18 ) = -27.0*xi*eta;
268
269 adNdXi( 0, 19 ) = -27.0*eta*zeta;
270 adNdXi( 1, 19 ) = 27.0*zeta*(1.0-(xi+2.0*eta+zeta));
271 adNdXi( 2, 19 ) = 27.0*eta*(1.0-(xi+eta+2.0*zeta));
272 }
273
274//------------------------------------------------------------------------------
275
276 template<>
277 void
280 const Vector< real > & aXi,
281 Matrix< real > & ad2NdXi2 ) const
282 {
283 const real xi = aXi( 0 );
284 const real eta = aXi( 1 );
285 const real zeta = aXi( 2 );
286
287 ad2NdXi2.set_size( 6, 20 );
288
289 ad2NdXi2( 0, 0 ) = 27.0*xi-9.0;
290 ad2NdXi2( 1, 0 ) = 0.0;
291 ad2NdXi2( 2, 0 ) = 0.0;
292 ad2NdXi2( 3, 0 ) = 0.0;
293 ad2NdXi2( 4, 0 ) = 0.0;
294 ad2NdXi2( 5, 0 ) = 0.0;
295
296 ad2NdXi2( 0, 1 ) = 0.0;
297 ad2NdXi2( 1, 1 ) = 0.0;
298 ad2NdXi2( 2, 1 ) = 27.0*zeta-9.0;
299 ad2NdXi2( 3, 1 ) = 0.0;
300 ad2NdXi2( 4, 1 ) = 0.0;
301 ad2NdXi2( 5, 1 ) = 0.0;
302
303 ad2NdXi2( 0, 2 ) = 0.0;
304 ad2NdXi2( 1, 2 ) = 27.0*eta-9.0;
305 ad2NdXi2( 2, 2 ) = 0.0;
306 ad2NdXi2( 3, 2 ) = 0.0;
307 ad2NdXi2( 4, 2 ) = 0.0;
308 ad2NdXi2( 5, 2 ) = 0.0;
309
310 ad2NdXi2( 0, 3 ) = 18.0-27.0*(xi+eta+zeta);
311 ad2NdXi2( 1, 3 ) = 18.0-27.0*(xi+eta+zeta);
312 ad2NdXi2( 2, 3 ) = 18.0-27.0*(xi+eta+zeta);
313 ad2NdXi2( 3, 3 ) = 18.0-27.0*(xi+eta+zeta);
314 ad2NdXi2( 4, 3 ) = 18.0-27.0*(xi+eta+zeta);
315 ad2NdXi2( 5, 3 ) = 18.0-27.0*(xi+eta+zeta);
316
317 ad2NdXi2( 0, 4 ) = 27.0*zeta;
318 ad2NdXi2( 1, 4 ) = 0.0;
319 ad2NdXi2( 2, 4 ) = 0.0;
320 ad2NdXi2( 3, 4 ) = 0.0;
321 ad2NdXi2( 4, 4 ) = 27.0*xi-4.5;
322 ad2NdXi2( 5, 4 ) = 0.0;
323
324 ad2NdXi2( 0, 5 ) = 0.0;
325 ad2NdXi2( 1, 5 ) = 0.0;
326 ad2NdXi2( 2, 5 ) = 27.0*xi;
327 ad2NdXi2( 3, 5 ) = 0.0;
328 ad2NdXi2( 4, 5 ) = 27.0*zeta-4.5;
329 ad2NdXi2( 5, 5 ) = 0.0;
330
331 ad2NdXi2( 0, 6 ) = 0.0;
332 ad2NdXi2( 1, 6 ) = 0.0;
333 ad2NdXi2( 2, 6 ) = 27.0*eta;
334 ad2NdXi2( 3, 6 ) = 27.0*zeta-4.5;
335 ad2NdXi2( 4, 6 ) = 0.0;
336 ad2NdXi2( 5, 6 ) = 0.0;
337
338 ad2NdXi2( 0, 7 ) = 0.0;
339 ad2NdXi2( 1, 7 ) = 27.0*zeta;
340 ad2NdXi2( 2, 7 ) = 0.0;
341 ad2NdXi2( 3, 7 ) = 27.0*eta-4.5;
342 ad2NdXi2( 4, 7 ) = 0.0;
343 ad2NdXi2( 5, 7 ) = 0.0;
344
345 ad2NdXi2( 0, 8 ) = 0.0;
346 ad2NdXi2( 1, 8 ) = 27.0*xi;
347 ad2NdXi2( 2, 8 ) = 0.0;
348 ad2NdXi2( 3, 8 ) = 0.0;
349 ad2NdXi2( 4, 8 ) = 0.0;
350 ad2NdXi2( 5, 8 ) = 27.0*eta-4.5;
351
352 ad2NdXi2( 0, 9 ) = 27.0*eta;
353 ad2NdXi2( 1, 9 ) = 0.0;
354 ad2NdXi2( 2, 9 ) = 0.0;
355 ad2NdXi2( 3, 9 ) = 0.0;
356 ad2NdXi2( 4, 9 ) = 0.0;
357 ad2NdXi2( 5, 9 ) = 27.0*xi-4.5;
358
359
360 ad2NdXi2( 0, 10 ) = 81.0*xi+54.0*(eta+zeta)-45.0;
361 ad2NdXi2( 1, 10 ) = 27.0*xi;
362 ad2NdXi2( 2, 10 ) = 27.0*xi;
363 ad2NdXi2( 3, 10 ) = 27.0*xi;
364 ad2NdXi2( 4, 10 ) = 54.0*xi+27.0*(eta+zeta)-22.5;
365 ad2NdXi2( 5, 10 ) = 54.0*xi+27.0*(eta+zeta)-22.5;
366
367 ad2NdXi2( 0, 11 ) = 36.0-81.0*xi-27.0*(eta+zeta);
368 ad2NdXi2( 1, 11 ) = 0.0;
369 ad2NdXi2( 2, 11 ) = 0.0;
370 ad2NdXi2( 3, 11 ) = 0.0;
371 ad2NdXi2( 4, 11 ) = 4.5-27.0*xi;
372 ad2NdXi2( 5, 11 ) = 4.5-27.0*xi;
373
374 ad2NdXi2( 0, 12 ) = 27.0*eta;
375 ad2NdXi2( 1, 12 ) = 54.0*(xi+zeta)+81.0*eta-45.0;
376 ad2NdXi2( 2, 12 ) = 27.0*eta;
377 ad2NdXi2( 3, 12 ) = 27.0*(xi+zeta)+54.0*eta-22.5;
378 ad2NdXi2( 4, 12 ) = 27.0*eta;
379 ad2NdXi2( 5, 12 ) = 27.0*(xi+zeta)+54.0*eta-22.5;
380
381 ad2NdXi2( 0, 13 ) = 0.0;
382 ad2NdXi2( 1, 13 ) = 36.0-27.0*(xi+zeta)-81.0*eta;
383 ad2NdXi2( 2, 13 ) = 0.0;
384 ad2NdXi2( 3, 13 ) = 4.5-27.0*eta;
385 ad2NdXi2( 4, 13 ) = 0.0;
386 ad2NdXi2( 5, 13 ) = 4.5-27.0*eta;
387
388 ad2NdXi2( 0, 14 ) = 27.0*zeta;
389 ad2NdXi2( 1, 14 ) = 27.0*zeta;
390 ad2NdXi2( 2, 14 ) = 54.0*(xi+eta)+81.0*zeta-45.0;
391 ad2NdXi2( 3, 14 ) = 27.0*(xi+eta)+54.0*zeta-22.5;
392 ad2NdXi2( 4, 14 ) = 27.0*(xi+eta)+54.0*zeta-22.5;
393 ad2NdXi2( 5, 14 ) = 27.0*zeta;
394
395 ad2NdXi2( 0, 15 ) = 0.0;
396 ad2NdXi2( 1, 15 ) = 0.0;
397 ad2NdXi2( 2, 15 ) = 36.0-27.0*(xi+eta)-81.0*zeta;
398 ad2NdXi2( 3, 15 ) = 4.5-27.0*zeta;
399 ad2NdXi2( 4, 15 ) = 4.5-27.0*zeta;
400 ad2NdXi2( 5, 15 ) = 0.0;
401
402 ad2NdXi2( 0, 16 ) = 0.0;
403 ad2NdXi2( 1, 16 ) = 0.0;
404 ad2NdXi2( 2, 16 ) = 0.0;
405 ad2NdXi2( 3, 16 ) = 27.0*xi;
406 ad2NdXi2( 4, 16 ) = 27.0*eta;
407 ad2NdXi2( 5, 16 ) = 27.0*zeta;
408
409 ad2NdXi2( 0, 17 ) = -54.0*zeta;
410 ad2NdXi2( 1, 17 ) = 0.0;
411 ad2NdXi2( 2, 17 ) = -54.0*xi;
412 ad2NdXi2( 3, 17 ) = -27.0*xi;
413 ad2NdXi2( 4, 17 ) = 27.0-54.0*(xi+zeta)-27.0*eta;
414 ad2NdXi2( 5, 17 ) = -27.0*zeta;
415
416 ad2NdXi2( 0, 18 ) = -54.0*eta;
417 ad2NdXi2( 1, 18 ) = -54.0*xi;
418 ad2NdXi2( 2, 18 ) = 0.0;
419 ad2NdXi2( 3, 18 ) = -27.0*xi;
420 ad2NdXi2( 4, 18 ) = -27.0*eta;
421 ad2NdXi2( 5, 18 ) = 27.0-54.0*(xi+eta)-27.0*zeta;
422
423 ad2NdXi2( 0, 19 ) = 0.0;
424 ad2NdXi2( 1, 19 ) = -54.0*zeta;
425 ad2NdXi2( 2, 19 ) = -54.0*eta;
426 ad2NdXi2( 3, 19 ) = 27.0-27.0*xi-54.0*(eta+zeta);
427 ad2NdXi2( 4, 19 ) = -27.0*eta;
428 ad2NdXi2( 5, 19 ) = -27.0*zeta;
429 }
430
431//------------------------------------------------------------------------------
432 }
433}
434#endif //BELFEM_CL_IF_TET20_HPP
void set_size(const size_t aNumRows, const size_t aNumCols)
Definition cl_AR_Matrix.hpp:186
shape function templated class G : Geometry T : Type D : Dimension B : Number of Basis
Definition cl_IF_InterpolationFunctionTemplate.hpp:25
void param_coords(Matrix< real > &aXiHat) const override
returns a matrix containing the parameter coordinates of the nodes < number of dimensions x number of...
Definition cl_IF_InterpolationFunctionTemplate.hpp:49
InterpolationOrder interpolation_order() const override
returns the interpolation order
Definition cl_IF_InterpolationFunctionTemplate.hpp:145
void d2NdXi2(const Vector< real > &aXi, Matrix< real > &ad2NdXi2) const override
calculates the second derivative of the shape function in parameter space
Definition cl_IF_InterpolationFunctionTemplate.hpp:110
void dNdXi(const Vector< real > &aXi, Matrix< real > &adNdXi) const override
calculates the first derivative of the shape function in parameter space
Definition cl_IF_InterpolationFunctionTemplate.hpp:89
void N(const Vector< real > &aXi, Matrix< real > &aN) const override
evaluates the shape function at a given point
Definition cl_IF_InterpolationFunctionTemplate.hpp:68
Definition cl_IFB_LINE3.hpp:21
USER GUIDES:
Definition cl_Capacitor.cpp:16
@ LAGRANGE
Definition Mesh_Enums.hpp:100
ElementType element_type(const std::string &aStr)
Definition Mesh_Enums.hpp:370
ElementType
Element types.
Definition Mesh_Enums.hpp:27
@ TET20
Definition Mesh_Enums.hpp:54
InterpolationOrder
Definition Mesh_Enums.hpp:85
@ CUBIC
Definition Mesh_Enums.hpp:90
double real
Definition typedefs.hpp:36
@ TET
Definition Mesh_Enums.hpp:75