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