BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
cl_IF_QUAD16.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_QUAD16_HPP
13#define BELFEM_CL_IF_QUAD16_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::QUAD, InterpolationType::LAGRANGE, 2, 16 >
31
32//------------------------------------------------------------------------------
33
34 template<>
35 ElementType
36 InterpolationFunctionTemplate<
37 GeometryType::QUAD, InterpolationType::LAGRANGE, 2, 16 >
38 ::element_type() const
39 {
40 return ElementType::QUAD16;
41 }
42
43//------------------------------------------------------------------------------
44
45 template<>
46 void
47 InterpolationFunctionTemplate<
48 GeometryType::QUAD, InterpolationType::LAGRANGE, 2, 16 >
49 ::param_coords( Matrix< real > & aXiHat ) const
50 {
51 aXiHat.set_size( 2, 16 );
52
53 const real c = 1.0/3.0;
54 aXiHat( 0, 0 ) = -1.000000;
55 aXiHat( 1, 0 ) = -1.000000;
56
57 aXiHat( 0, 1 ) = 1.000000;
58 aXiHat( 1, 1 ) = -1.000000;
59
60 aXiHat( 0, 2 ) = 1.000000;
61 aXiHat( 1, 2 ) = 1.000000;
62
63 aXiHat( 0, 3 ) = -1.000000;
64 aXiHat( 1, 3 ) = 1.000000;
65
66 aXiHat( 0, 4 ) = -c;
67 aXiHat( 1, 4 ) = -1.000000;
68
69 aXiHat( 0, 5 ) = c;
70 aXiHat( 1, 5 ) = -1.000000;
71
72 aXiHat( 0, 6 ) = 1.000000;
73 aXiHat( 1, 6 ) = -c;
74
75 aXiHat( 0, 7 ) = 1.000000;
76 aXiHat( 1, 7 ) = c;
77
78 aXiHat( 0, 8 ) = c;
79 aXiHat( 1, 8 ) = 1.000000;
80
81 aXiHat( 0, 9 ) = -c;
82 aXiHat( 1, 9 ) = 1.000000;
83
84 aXiHat( 0, 10 ) = -1.000000;
85 aXiHat( 1, 10 ) = c;
86
87 aXiHat( 0, 11 ) = -1.000000;
88 aXiHat( 1, 11 ) = -c;
89
90 aXiHat( 0, 12 ) = -c;
91 aXiHat( 1, 12 ) = -c;
92
93 aXiHat( 0, 13 ) = c;
94 aXiHat( 1, 13 ) = -c;
95
96 aXiHat( 0, 14 ) = c;
97 aXiHat( 1, 14 ) = c;
98
99 aXiHat( 0, 15 ) = -c;
100 aXiHat( 1, 15 ) = c;
101 }
102
103//------------------------------------------------------------------------------
104
105 template<>
106 void
107 InterpolationFunctionTemplate<
108 GeometryType::QUAD, InterpolationType::LAGRANGE, 2, 16 >::N(
109 const Vector< real > & aXi,
110 Matrix< real > & aN ) const
111 {
112 const real xi = aXi( 0 );
113 const real eta = aXi( 1 );
114
115 const real a0 = ( xi*( 1.0 + 9.0 * xi * ( 1.0 - xi ) ) - 1.0 )*0.0625;
116 const real a1 = ( 9.0 - xi * ( 27.0 + xi*( 9.0 - 27.0*xi ) ) )*0.0625;
117 const real a2 = ( 9.0 + xi * ( 27.0 - xi*( 9.0 + 27.0*xi ) ) )*0.0625;
118 const real a3 = ( -xi*( 1.0 - 9.0 * xi * ( 1.0 + xi ) ) - 1.0 )*0.0625;
119
120 const real b0 = ( eta*( 1.0 + 9.0 * eta * ( 1.0 - eta ) ) - 1.0 )*0.0625;
121 const real b1 = ( 9.0 - eta * ( 27.0 + eta*( 9.0 - 27.0*eta ) ) )*0.0625;
122 const real b2 = ( 9.0 + eta * ( 27.0 - eta*( 9.0 + 27.0*eta ) ) )*0.0625;
123 const real b3 = ( -eta*( 1.0 - 9.0 * eta * ( 1.0 + eta ) ) - 1.0 )*0.0625;
124
125 // populate matrix with values
126 aN.set_size(1,16);
127 aN( 0, 0 ) = a0*b0;
128 aN( 0, 1 ) = a3*b0;
129 aN( 0, 2 ) = a3*b3;
130 aN( 0, 3 ) = a0*b3;
131 aN( 0, 4 ) = a1*b0;
132 aN( 0, 5 ) = a2*b0;
133 aN( 0, 6 ) = a3*b1;
134 aN( 0, 7 ) = a3*b2;
135 aN( 0, 8 ) = a2*b3;
136 aN( 0, 9 ) = a1*b3;
137 aN( 0, 10 ) = a0*b2;
138 aN( 0, 11 ) = a0*b1;
139 aN( 0, 12 ) = a1*b1;
140 aN( 0, 13 ) = a2*b1;
141 aN( 0, 14 ) = a2*b2;
142 aN( 0, 15 ) = a1*b2;
143 }
144
145//------------------------------------------------------------------------------
146
147 template<>
148 void
149 InterpolationFunctionTemplate<
150 GeometryType::QUAD, InterpolationType::LAGRANGE, 2, 16 >::dNdXi(
151 const Vector< real > & aXi,
152 Matrix< real > & adNdXi ) const
153 {
154 const real xi = aXi( 0 );
155 const real eta = aXi( 1 );
156
157 // often used parameters
158 const real a0 = ( xi*( 1.0 + 9.0 * xi * ( 1.0 - xi ) ) - 1.0 ) * 0.0625;
159 const real a1 = ( 9.0 - xi * ( 27.0 + xi*( 9.0 - 27.0*xi ) ) ) * 0.0625;
160 const real a2 = ( 9.0 + xi * ( 27.0 - xi*( 9.0 + 27.0*xi ) ) ) * 0.0625;
161 const real a3 = ( -xi*( 1.0 - 9.0 * xi * ( 1.0 + xi ) ) - 1.0 ) * 0.0625;
162
163 const real b0 = ( eta*( 1.0 + 9.0 * eta * ( 1.0 - eta ) ) - 1.0 ) * 0.0625;
164 const real b1 = ( 9.0 - eta * ( 27.0 + eta*( 9.0 - 27.0*eta ) ) ) * 0.0625;
165 const real b2 = ( 9.0 + eta * ( 27.0 - eta*( 9.0 + 27.0*eta ) ) ) * 0.0625;
166 const real b3 = ( -eta*( 1.0 - 9.0 * eta * ( 1.0 + eta ) ) - 1.0 ) * 0.0625;
167
168 const real da0 = ( 1.0 + xi*( 18.0 - 27.0*xi )) * 0.0625;
169 const real da1 = ( -27.0 - xi*( 18.0 - 81.0*xi )) * 0.0625;
170 const real da2 = ( 27.0 - xi*( 18.0 + 81.0*xi )) * 0.0625;
171 const real da3 = ( -1.0 + xi*( 18.0 + 27.0*xi )) * 0.0625;
172
173 const real db0 = ( 1.0 + eta*( 18.0 - 27.0*eta )) * 0.0625;
174 const real db1 = ( -27.0 - eta*( 18.0 - 81.0*eta )) * 0.0625;
175 const real db2 = ( 27.0 - eta*( 18.0 + 81.0*eta )) * 0.0625;
176 const real db3 = ( -1.0 + eta*( 18.0 + 27.0*eta )) * 0.0625;
177
178 // populate output matrix
179 adNdXi.set_size( 2, 16 );
180
181 adNdXi( 0, 0 ) = da0*b0;
182 adNdXi( 1, 0 ) = a0*db0;
183
184 adNdXi( 0, 1 ) = da3*b0;
185 adNdXi( 1, 1 ) = a3*db0;
186
187 adNdXi( 0, 2 ) = da3*b3;
188 adNdXi( 1, 2 ) = a3*db3;
189
190 adNdXi( 0, 3 ) = da0*b3;
191 adNdXi( 1, 3 ) = a0*db3;
192
193 adNdXi( 0, 4 ) = da1*b0;
194 adNdXi( 1, 4 ) = a1*db0;
195
196 adNdXi( 0, 5 ) = da2*b0;
197 adNdXi( 1, 5 ) = a2*db0;
198
199 adNdXi( 0, 6 ) = da3*b1;
200 adNdXi( 1, 6 ) = a3*db1;
201
202 adNdXi( 0, 7 ) = da3*b2;
203 adNdXi( 1, 7 ) = a3*db2;
204
205 adNdXi( 0, 8 ) = da2*b3;
206 adNdXi( 1, 8 ) = a2*db3;
207
208 adNdXi( 0, 9 ) = da1*b3;
209 adNdXi( 1, 9 ) = a1*db3;
210
211 adNdXi( 0, 10 ) = da0*b2;
212 adNdXi( 1, 10 ) = a0*db2;
213
214 adNdXi( 0, 11 ) = da0*b1;
215 adNdXi( 1, 11 ) = a0*db1;
216
217 adNdXi( 0, 12 ) = da1*b1;
218 adNdXi( 1, 12 ) = a1*db1;
219
220 adNdXi( 0, 13 ) = da2*b1;
221 adNdXi( 1, 13 ) = a2*db1;
222
223 adNdXi( 0, 14 ) = da2*b2;
224 adNdXi( 1, 14 ) = a2*db2;
225
226 adNdXi( 0, 15 ) = da1*b2;
227 adNdXi( 1, 15 ) = a1*db2;
228 }
229
230//------------------------------------------------------------------------------
231
232 template<>
233 void
234 InterpolationFunctionTemplate<
235 GeometryType::QUAD, InterpolationType::LAGRANGE, 2, 16 >::d2NdXi2(
236 const Vector< real > & aXi,
237 Matrix< real > & ad2NdXi2 ) const
238 {
239 // unpack xi and eta from input vector
240 const real xi = aXi( 0 );
241 const real eta = aXi( 1 );
242
243 // often used parameters
244 const real a0 = ( xi*( 1.0 + 9.0 * xi * ( 1.0 - xi ) ) - 1.0 ) * 0.0625;
245 const real a1 = ( 9.0 - xi * ( 27.0 + xi*( 9.0 - 27.0*xi ) ) ) * 0.0625;
246 const real a2 = ( 9.0 + xi * ( 27.0 - xi*( 9.0 + 27.0*xi ) ) ) * 0.0625;
247 const real a3 = ( -xi*( 1.0 - 9.0 * xi * ( 1.0 + xi ) ) - 1.0 ) * 0.0625;
248
249 const real b0 = ( eta*( 1.0 + 9.0 * eta * ( 1.0 - eta ) ) - 1.0 ) * 0.0625;
250 const real b1 = ( 9.0 - eta * ( 27.0 + eta*( 9.0 - 27.0*eta ) ) ) * 0.0625;
251 const real b2 = ( 9.0 + eta * ( 27.0 - eta*( 9.0 + 27.0*eta ) ) ) * 0.0625;
252 const real b3 = ( -eta*( 1.0 - 9.0 * eta * ( 1.0 + eta ) ) - 1.0 ) * 0.0625;
253
254 const real da0 = ( 1.0 + xi*( 18.0 - 27.0*xi )) * 0.0625;
255 const real da1 = ( -27.0 - xi*( 18.0 - 81.0*xi )) * 0.0625;
256 const real da2 = ( 27.0 - xi*( 18.0 + 81.0*xi )) * 0.0625;
257 const real da3 = ( -1.0 + xi*( 18.0 + 27.0*xi )) * 0.0625;
258
259 const real dda0 = ( 18.0 - 54.0*xi ) * 0.0625;
260 const real dda1 = ( 162.0*xi - 18.0 ) * 0.0625;
261 const real dda2 = ( - 162.0*xi - 18.0 ) * 0.0625;
262 const real dda3 = ( 54.0*xi + 18.0 ) * 0.0625;
263
264 const real db0 = ( 1.0 + eta*( 18.0 - 27.0*eta )) * 0.0625;
265 const real db1 = ( -27.0 - eta*( 18.0 - 81.0*eta )) * 0.0625;
266 const real db2 = ( 27.0 - eta*( 18.0 + 81.0*eta )) * 0.0625;
267 const real db3 = ( -1.0 + eta*( 18.0 + 27.0*eta )) * 0.0625;
268
269 const real ddb0 = ( 18.0 - 54.0*eta ) * 0.0625;
270 const real ddb1 = ( 162.0*eta - 18.0 ) * 0.0625;
271 const real ddb2 = ( - 162.0*eta - 18.0 ) * 0.0625;
272 const real ddb3 = ( 54.0*eta + 18.0 ) * 0.0625;
273
274 ad2NdXi2.set_size( 3, 16 );
275 ad2NdXi2( 0, 0 ) = dda0*b0;
276 ad2NdXi2( 1, 0 ) = a0*ddb0;
277 ad2NdXi2( 2, 0 ) = da0*db0;
278
279 ad2NdXi2( 0, 1 ) = dda3*b0;
280 ad2NdXi2( 1, 1 ) = a3*ddb0;
281 ad2NdXi2( 2, 1 ) = da3*db0;
282
283 ad2NdXi2( 0, 2 ) = dda3*b3;
284 ad2NdXi2( 1, 2 ) = a3*ddb3;
285 ad2NdXi2( 2, 2 ) = da3*db3;
286
287 ad2NdXi2( 0, 3 ) = dda0*b3;
288 ad2NdXi2( 1, 3 ) = a0*ddb3;
289 ad2NdXi2( 2, 3 ) = da0*db3;
290
291 ad2NdXi2( 0, 4 ) = dda1*b0;
292 ad2NdXi2( 1, 4 ) = a1*ddb0;
293 ad2NdXi2( 2, 4 ) = da1*db0;
294
295 ad2NdXi2( 0, 5 ) = dda2*b0;
296 ad2NdXi2( 1, 5 ) = a2*ddb0;
297 ad2NdXi2( 2, 5 ) = da2*db0;
298
299 ad2NdXi2( 0, 6 ) = dda3*b1;
300 ad2NdXi2( 1, 6 ) = a3*ddb1;
301 ad2NdXi2( 2, 6 ) = da3*db1;
302
303 ad2NdXi2( 0, 7 ) = dda3*b2;
304 ad2NdXi2( 1, 7 ) = a3*ddb2;
305 ad2NdXi2( 2, 7 ) = da3*db2;
306
307 ad2NdXi2( 0, 8 ) = dda2*b3;
308 ad2NdXi2( 1, 8 ) = a2*ddb3;
309 ad2NdXi2( 2, 8 ) = da2*db3;
310
311 ad2NdXi2( 0, 9 ) = dda1*b3;
312 ad2NdXi2( 1, 9 ) = a1*ddb3;
313 ad2NdXi2( 2, 9 ) = da1*db3;
314
315 ad2NdXi2( 0, 10 ) = dda0*b2;
316 ad2NdXi2( 1, 10 ) = a0*ddb2;
317 ad2NdXi2( 2, 10 ) = da0*db2;
318
319 ad2NdXi2( 0, 11 ) = dda0*b1;
320 ad2NdXi2( 1, 11 ) = a0*ddb1;
321 ad2NdXi2( 2, 11 ) = da0*db1;
322
323 ad2NdXi2( 0, 12 ) = dda1*b1;
324 ad2NdXi2( 1, 12 ) = a1*ddb1;
325 ad2NdXi2( 2, 12 ) = da1*db1;
326
327 ad2NdXi2( 0, 13 ) = dda2*b1;
328 ad2NdXi2( 1, 13 ) = a2*ddb1;
329 ad2NdXi2( 2, 13 ) = da2*db1;
330
331 ad2NdXi2( 0, 14 ) = dda2*b2;
332 ad2NdXi2( 1, 14 ) = a2*ddb2;
333 ad2NdXi2( 2, 14 ) = da2*db2;
334
335 ad2NdXi2( 0, 15 ) = dda1*b2;
336 ad2NdXi2( 1, 15 ) = a1*ddb2;
337 ad2NdXi2( 2, 15 ) = da1*db2;
338 }
339
340//------------------------------------------------------------------------------
341 }
342}
343
344#endif //BELFEM_CL_IF_QUAD16_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::QUAD16
@ QUAD16
Definition Mesh_Enums.hpp:56
belfem::InterpolationOrder
InterpolationOrder
Definition Mesh_Enums.hpp:85
belfem::InterpolationOrder::CUBIC
@ CUBIC
Definition Mesh_Enums.hpp:90
belfem::real
double real
Definition typedefs.hpp:36
belfem::GeometryType::QUAD
@ QUAD
Definition Mesh_Enums.hpp:74