BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
cl_IF_PYRA13.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_PYRA13_HPP
13#define BELFEM_CL_IF_PYRA13_HPP
14#include "cl_IF_InterpolationFunctionTemplate.hpp"
15
16namespace belfem
17{
18 namespace fem
19 {
20//------------------------------------------------------------------------------
21
22 template<>
23 InterpolationOrder
24 InterpolationFunctionTemplate<
25 GeometryType::PYRA, InterpolationType::LAGRANGE, 3, 13 >
30
31//------------------------------------------------------------------------------
32
33 template<>
34 ElementType
35 InterpolationFunctionTemplate<
36 GeometryType::PYRA, InterpolationType::LAGRANGE, 3, 13 >
37 ::element_type() const
38 {
39 return ElementType::PYRA13;
40 }
41
42//------------------------------------------------------------------------------
43
44 template<>
45 void
46 InterpolationFunctionTemplate<
47 GeometryType::PYRA, InterpolationType::LAGRANGE, 3, 13 >
48 ::param_coords( Matrix< real > & aXiHat ) const
49 {
50 aXiHat.set_size( 3, 13 );
51
52 aXiHat( 0, 0 ) = -1.0;
53 aXiHat( 1, 0 ) = -1.0;
54 aXiHat( 2, 0 ) = 0.0;
55
56 aXiHat( 0, 1 ) = 1.0;
57 aXiHat( 1, 1 ) = -1.0;
58 aXiHat( 2, 1 ) = 0.0;
59
60 aXiHat( 0, 2 ) = 1.0;
61 aXiHat( 1, 2 ) = 1.0;
62 aXiHat( 2, 2 ) = 0.0;
63
64 aXiHat( 0, 3 ) = -1.0;
65 aXiHat( 1, 3 ) = 1.0;
66 aXiHat( 2, 3 ) = 0.0;
67
68 aXiHat( 0, 4 ) = 0.0;
69 aXiHat( 1, 4 ) = 0.0;
70 aXiHat( 2, 4 ) = 1.0;
71
72 aXiHat( 0, 5 ) = 0.0 ;
73 aXiHat( 1, 5 ) = -1.0 ;
74 aXiHat( 2, 5 ) = 0.0 ;
75
76 aXiHat( 0, 6 ) = 1.0 ;
77 aXiHat( 1, 6 ) = 0.0 ;
78 aXiHat( 2, 6 ) = 0.0 ;
79
80 aXiHat( 0, 7 ) = 0.0 ;
81 aXiHat( 1, 7 ) = 1.0 ;
82 aXiHat( 2, 7 ) = 0.0 ;
83
84 aXiHat( 0, 8 ) = -1.0 ;
85 aXiHat( 1, 8 ) = 0.0 ;
86 aXiHat( 2, 8 ) = 0.0 ;
87
88 aXiHat( 0, 9 ) = -0.5 ;
89 aXiHat( 1, 9 ) = -0.5 ;
90 aXiHat( 2, 9 ) = 0.5 ;
91
92 aXiHat( 0, 10 ) = 0.5 ;
93 aXiHat( 1, 10 ) = -0.5 ;
94 aXiHat( 2, 10 ) = 0.5 ;
95
96 aXiHat( 0, 11 ) = 0.5 ;
97 aXiHat( 1, 11 ) = 0.5 ;
98 aXiHat( 2, 11 ) = 0.5 ;
99
100 aXiHat( 0, 12 ) = -0.5 ;
101 aXiHat( 1, 12 ) = 0.5 ;
102 aXiHat( 2, 12 ) = 0.5 ;
103
104 }
105
106//------------------------------------------------------------------------------
107
108 template<>
109 void
110 InterpolationFunctionTemplate<
111 GeometryType::PYRA, InterpolationType::LAGRANGE, 3, 13 >
112 ::N(
113 const Vector< real > & aXi,
114 Matrix< real > & aN ) const
115 {
116 const real xi = aXi( 0 );
117 const real eta = aXi( 1 );
118 const real zeta = aXi( 2 );
119
120 real xi2 = xi * xi;
121 real eta2 = eta * eta;
122 real zeta2 = zeta * zeta;
123 real phi = eta * zeta;
124 real psi = xi * zeta;
125 real chi = xi * eta;
126 real alpha = eta * xi2;
127 real beta = xi * eta2;
128 real gamma = chi * zeta;
129 real gamma2 = gamma + gamma ;
130
131 aN.set_size( 1, 13 );
132
133 aN( 0, 0 ) = 0.25 * ( zeta + chi - alpha - beta + xi2 + eta2 ) + 0.125 * ( phi + psi ) - 0.5*gamma - 0.25 ;
134 aN( 0, 1 ) = 0.25 * ( zeta - chi - alpha + beta + xi2 + eta2 ) + 0.125 * ( phi - psi ) + 0.5*gamma - 0.25 ;
135 aN( 0, 2 ) = 0.25 * ( zeta + chi + alpha + beta + xi2 + eta2 ) - 0.125 * ( phi + psi ) - 0.5*gamma - 0.25 ;
136 aN( 0, 3 ) = 0.25 * ( zeta - chi + alpha - beta + xi2 + eta2 ) - 0.125 * ( phi - psi ) + 0.5*gamma - 0.25 ;
137
138 aN( 0, 4 ) = zeta2 + zeta2 - zeta ;
139
140 aN( 0, 5 ) = 0.5*( zeta2 - xi2 + alpha - eta ) + 0.75 * phi - zeta + 0.5 ;
141 aN( 0, 6 ) = 0.5*( zeta2 - eta2 - beta + xi ) - 0.75 * psi - zeta + 0.5 ;
142 aN( 0, 7 ) = 0.5*( zeta2 - xi2 - alpha + eta ) - 0.75 * phi - zeta + 0.5 ;
143 aN( 0, 8 ) = 0.5*( zeta2 - eta2 + beta - xi ) + 0.75 * psi - zeta + 0.5 ;
144
145 aN( 0, 9 ) = zeta - phi - psi - zeta2 + gamma2 ;
146 aN( 0, 10 ) = zeta - phi + psi - zeta2 - gamma2 ;
147 aN( 0, 11 ) = zeta + phi + psi - zeta2 + gamma2 ;
148 aN( 0, 12 ) = zeta + phi - psi - zeta2 - gamma2 ;
149 }
150
151//------------------------------------------------------------------------------
152
153 template<>
154 void
155 InterpolationFunctionTemplate<
156 GeometryType::PYRA, InterpolationType::LAGRANGE, 3, 13 >
157 ::dNdXi(
158 const Vector< real > & aXi,
159 Matrix< real > & adNdXi ) const
160 {
161 const real xi = aXi( 0 );
162 const real eta = aXi( 1 );
163 const real zeta = aXi( 2 );
164
165 real xi2 = xi * xi;
166 real eta2 = eta * eta;
167
168 real phi = eta * zeta;
169 real psi = xi * zeta;
170 real chi = xi * eta;
171
172 real phi2 = phi + phi ;
173 real psi2 = psi + psi ;
174 real zeta8 = 0.125 * zeta ;
175
176 adNdXi.set_size( 3, 13 );
177
178 adNdXi( 0, 0 ) = 0.5 * ( xi - chi - phi ) + 0.25 * ( eta - eta2 ) + zeta8 ;
179 adNdXi( 1, 0 ) = 0.5 * ( eta - chi - psi ) + 0.25 * ( xi - xi2 ) + zeta8 ;
180 adNdXi( 2, 0 ) = 0.25 + 0.125 * ( xi + eta ) - 0.5 * chi ;
181
182 adNdXi( 0, 1 ) = 0.5 * ( xi - chi + phi ) + 0.25 * ( eta2 - eta ) - zeta8 ;
183 adNdXi( 1, 1 ) = 0.5 * ( eta + chi + psi ) - 0.25 * ( xi + xi2 ) + zeta8 ;
184 adNdXi( 2, 1 ) = 0.25 + 0.125 * ( eta - xi ) + 0.5 * chi ;
185
186 adNdXi( 0, 2 ) = 0.5 * ( xi + chi - phi ) + 0.25 * ( eta2 + eta ) - zeta8 ;
187 adNdXi( 1, 2 ) = 0.5 * ( eta + chi - psi ) + 0.25 * ( xi + xi2 ) - zeta8 ;
188 adNdXi( 2, 2 ) = 0.25 - 0.125 * ( xi + eta ) - 0.5 * chi ;
189
190 adNdXi( 0, 3 ) = 0.5 * ( xi + chi + phi ) - 0.25 * ( eta2 + eta ) + zeta8 ;
191 adNdXi( 1, 3 ) = 0.5 * ( eta - chi + psi ) + 0.25 * ( xi2 - xi ) - zeta8 ;
192 adNdXi( 2, 3 ) = 0.25 + 0.125 * ( xi - eta ) + 0.5 * chi ;
193
194 adNdXi( 0, 4 ) = 0.0 ;
195 adNdXi( 1, 4 ) = 0.0 ;
196 adNdXi( 2, 4 ) = 4.0 * zeta - 1.0 ;
197
198 adNdXi( 0, 5 ) = chi - xi ;
199 adNdXi( 1, 5 ) = 0.5 * xi2 + 0.75 * zeta - 0.5 ;
200 adNdXi( 2, 5 ) = zeta + 0.75 * eta - 1.0 ;
201
202 adNdXi( 0, 6 ) = 0.5 - 0.5 * eta2 - 0.75 * zeta ;
203 adNdXi( 1, 6 ) = -eta - chi ;
204 adNdXi( 2, 6 ) = zeta - 0.75 * xi - 1.0 ;
205
206 adNdXi( 0, 7 ) = -chi - xi ;
207 adNdXi( 1, 7 ) = 0.5 - 0.5 * xi2 - 0.75 * zeta ;
208 adNdXi( 2, 7 ) = zeta - 0.75 * eta - 1.0 ;
209
210 adNdXi( 0, 8 ) = 0.5 * eta2 + 0.75 * zeta - 0.5 ;
211 adNdXi( 1, 8 ) = chi - eta ;
212 adNdXi( 2, 8 ) = zeta + 0.75 * xi - 1.0 ;
213
214 adNdXi( 0, 9 ) = phi2 - zeta ;
215 adNdXi( 1, 9 ) = psi2 - zeta ;
216 adNdXi( 2, 9 ) = 1.0 - xi - eta + 2.0 * (chi - zeta);
217
218 adNdXi( 0, 10 ) = zeta - phi2 ;
219 adNdXi( 1, 10 ) = -psi2 - zeta ;
220 adNdXi( 2, 10 ) = 1.0 + xi - eta - 2.0 * (chi + zeta);
221
222 adNdXi( 0, 11 ) = zeta + phi2 ;
223 adNdXi( 1, 11 ) = zeta + psi2 ;
224 adNdXi( 2, 11 ) = 1.0 + xi + eta + 2.0 * (chi - zeta);
225
226 adNdXi( 0, 12 ) = -phi2 - zeta ;
227 adNdXi( 1, 12 ) = zeta - psi2 ;
228 adNdXi( 2, 12 ) = 1.0 - xi + eta - 2.0 * (chi + zeta);
229 }
230
231//------------------------------------------------------------------------------
232
233 template<>
234 void
235 InterpolationFunctionTemplate<
236 GeometryType::PYRA, InterpolationType::LAGRANGE, 3, 13 >::d2NdXi2(
237 const Vector< real > & aXi,
238 Matrix< real > & ad2NdXi2 ) const
239 {
240
241 const real xi = aXi( 0 );
242 const real eta = aXi( 1 );
243 const real zeta = aXi( 2 );
244
245 const real hxi = 0.5 * xi ;
246 const real heta = 0.5 * eta ;
247 const real hzeta = 0.5 * zeta ;
248
249 const real xi2 = xi + xi ;
250 const real eta2 = eta + eta ;
251 const real zeta2 = zeta + zeta ;
252
253 ad2NdXi2.set_size( 6, 13, 0.0 );
254
255 ad2NdXi2( 0, 0 ) = 0.5 - heta ;
256 ad2NdXi2( 1, 0 ) = 0.5 - hxi ;
257 ad2NdXi2( 2, 0 ) = 0.0 ;
258 ad2NdXi2( 3, 0 ) = 0.125 - hxi ;
259 ad2NdXi2( 4, 0 ) = 0.125 - heta ;
260 ad2NdXi2( 5, 0 ) = 0.25 - hxi - hzeta - heta ;
261
262 ad2NdXi2( 0, 1 ) = 0.5 - heta ;
263 ad2NdXi2( 1, 1 ) = hxi + 0.5 ;
264 ad2NdXi2( 2, 1 ) = 0.0 ;
265 ad2NdXi2( 3, 1 ) = hxi + 0.125 ;
266 ad2NdXi2( 4, 1 ) = heta - 0.125 ;
267 ad2NdXi2( 5, 1 ) = heta - hxi + hzeta - 0.25 ;
268
269 ad2NdXi2( 0, 2 ) = heta + 0.5 ;
270 ad2NdXi2( 1, 2 ) = hxi + 0.5 ;
271 ad2NdXi2( 2, 2 ) = 0.0 ;
272 ad2NdXi2( 3, 2 ) = - hxi - 0.125 ;
273 ad2NdXi2( 4, 2 ) = - heta - 0.125 ;
274 ad2NdXi2( 5, 2 ) = heta + hxi - hzeta + 0.25 ;
275
276 ad2NdXi2( 0, 3 ) = heta + 0.5 ;
277 ad2NdXi2( 1, 3 ) = 0.5 - hxi ;
278 ad2NdXi2( 2, 3 ) = 0.0 ;
279 ad2NdXi2( 3, 3 ) = hxi - 0.125 ;
280 ad2NdXi2( 4, 3 ) = heta + 0.125 ;
281 ad2NdXi2( 5, 3 ) = hxi - heta + hzeta - 0.25 ;
282
283 ad2NdXi2( 0, 4 ) = 0.0 ;
284 ad2NdXi2( 1, 4 ) = 0.0 ;
285 ad2NdXi2( 2, 4 ) = 4.0 ;
286 ad2NdXi2( 3, 4 ) = 0.0 ;
287 ad2NdXi2( 4, 4 ) = 0.0 ;
288 ad2NdXi2( 5, 4 ) = 0.0 ;
289
290 ad2NdXi2( 0, 5 ) = eta - 1.0 ;
291 ad2NdXi2( 1, 5 ) = 0.0 ;
292 ad2NdXi2( 2, 5 ) = 1.0 ;
293 ad2NdXi2( 3, 5 ) = 0.75 ;
294 ad2NdXi2( 4, 5 ) = 0.0 ;
295 ad2NdXi2( 5, 5 ) = xi ;
296
297 ad2NdXi2( 0, 6 ) = 0.0 ;
298 ad2NdXi2( 1, 6 ) = - 1.0 - xi ;
299 ad2NdXi2( 2, 6 ) = 1.0 ;
300 ad2NdXi2( 3, 6 ) = 0.0 ;
301 ad2NdXi2( 4, 6 ) = -0.75 ;
302 ad2NdXi2( 5, 6 ) = -eta ;
303
304 ad2NdXi2( 0, 7 ) = - 1.0 - eta ;
305 ad2NdXi2( 1, 7 ) = 0 ;
306 ad2NdXi2( 2, 7 ) = 1.0 ;
307 ad2NdXi2( 3, 7 ) = -0.75 ;
308 ad2NdXi2( 4, 7 ) = 0.0 ;
309 ad2NdXi2( 5, 7 ) = - xi ;
310
311 ad2NdXi2( 0, 8 ) = 0.0 ;
312 ad2NdXi2( 1, 8 ) = xi - 1.0 ;
313 ad2NdXi2( 2, 8 ) = 1.0 ;
314 ad2NdXi2( 3, 8 ) = 0.0 ;
315 ad2NdXi2( 4, 8 ) = 0.75 ;
316 ad2NdXi2( 5, 8 ) = eta ;
317
318 ad2NdXi2( 0, 9 ) = 0.0 ;
319 ad2NdXi2( 1, 9 ) = 0.0 ;
320 ad2NdXi2( 2, 9 ) = -2.0 ;
321 ad2NdXi2( 3, 9 ) = xi2 - 1.0 ;
322 ad2NdXi2( 4, 9 ) = eta2 - 1.0 ;
323 ad2NdXi2( 5, 9 ) = zeta2 ;
324
325 ad2NdXi2( 0, 10 ) = 0.0 ;
326 ad2NdXi2( 1, 10 ) = 0.0 ;
327 ad2NdXi2( 2, 10 ) = -2.0 ;
328 ad2NdXi2( 3, 10 ) = -xi2 - 1.0 ;
329 ad2NdXi2( 4, 10 ) = 1.0 - eta2 ;
330 ad2NdXi2( 5, 10 ) = - zeta2 ;
331
332 ad2NdXi2( 0, 11 ) = 0.0 ;
333 ad2NdXi2( 1, 11 ) = 0.0 ;
334 ad2NdXi2( 2, 11 ) = -2.0 ;
335 ad2NdXi2( 3, 11 ) = xi2 + 1.0 ;
336 ad2NdXi2( 4, 11 ) = eta2 + 1.0 ;
337 ad2NdXi2( 5, 11 ) = zeta2 ;
338
339 ad2NdXi2( 0, 12 ) = 0.0 ;
340 ad2NdXi2( 1, 12 ) = 0.0 ;
341 ad2NdXi2( 2, 12 ) = -2.0 ;
342 ad2NdXi2( 3, 12 ) = 1.0 - xi2 ;
343 ad2NdXi2( 4, 12 ) = -1.0 - eta2 ;
344 ad2NdXi2( 5, 12 ) = - zeta2 ;
345
346 }
347
348//------------------------------------------------------------------------------
349 }
350}
351
352#endif //BELFEM_CL_IF_PYRA13_HPP
cl_IF_InterpolationFunctionTemplate.hpp
belfem::Matrix::set_size
void set_size(const size_t aNumRows, const size_t aNumCols)
Definition cl_AR_Matrix.hpp:186
belfem::fem::InterpolationFunctionTemplate
shape function templated class G : Geometry T : Type D : Dimension B : Number of Basis
Definition cl_IF_InterpolationFunctionTemplate.hpp:25
belfem::fem::InterpolationFunctionTemplate::param_coords
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
belfem::fem::InterpolationFunctionTemplate::interpolation_order
InterpolationOrder interpolation_order() const override
returns the interpolation order
Definition cl_IF_InterpolationFunctionTemplate.hpp:145
belfem::fem::InterpolationFunctionTemplate::d2NdXi2
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
belfem::fem::InterpolationFunctionTemplate::dNdXi
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
belfem::fem::InterpolationFunctionTemplate::N
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
belfem::fem
Definition cl_IFB_LINE3.hpp:21
belfem
USER GUIDES:
Definition cl_Capacitor.cpp:16
belfem::InterpolationType::LAGRANGE
@ LAGRANGE
Definition Mesh_Enums.hpp:100
belfem::element_type
ElementType element_type(const std::string &aStr)
Definition Mesh_Enums.hpp:370
belfem::ElementType
ElementType
Element types.
Definition Mesh_Enums.hpp:27
belfem::ElementType::PYRA13
@ PYRA13
Definition Mesh_Enums.hpp:47
belfem::MaterialProperty::gamma
@ gamma
Definition cl_Material.hpp:174
belfem::MaterialProperty::alpha
@ alpha
Definition cl_Material.hpp:161
belfem::InterpolationOrder
InterpolationOrder
Definition Mesh_Enums.hpp:85
belfem::InterpolationOrder::SERENDIPITY
@ SERENDIPITY
Definition Mesh_Enums.hpp:89
belfem::real
double real
Definition typedefs.hpp:36
belfem::GeometryType::PYRA
@ PYRA
Definition Mesh_Enums.hpp:78