BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
cl_IFG_TET10B.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_IFG_TET10B_HPP
13#define BELFEM_CL_IFG_TET10B_HPP
14
16
17namespace belfem
18{
19 namespace fem
20 {
21//------------------------------------------------------------------------------
22
30
31//------------------------------------------------------------------------------
32
33 template<>
41
42//------------------------------------------------------------------------------
43
44 template<>
49 {
50 return ElementType::TET10;
51 }
52
53//------------------------------------------------------------------------------
54
55
56 template<>
57 void
60 ::param_coords( Matrix <real> & aXiHat ) const
61 {
62 aXiHat.set_size( 3, 3 );
63
64 // point 0
65 aXiHat( 0, 0 ) = 0.0 ;
66 aXiHat( 0, 1 ) = 2./3. ;
67 aXiHat( 0, 2 ) = 1./6. ;
68
69 // point 1
70 aXiHat( 1, 0 ) = 0.0 ;
71 aXiHat( 1, 1 ) = 1./6. ;
72 aXiHat( 1, 2 ) = 2./3. ;
73
74 // point 2
75 aXiHat( 2, 0 ) = 0.0 ;
76 aXiHat( 2, 1 ) = 1./6. ;
77 aXiHat( 2, 2 ) = 1./6. ;
78
79 }
80
81//------------------------------------------------------------------------------
82
83 template<>
84 void
87 const Vector <real> & aXi,
88 Matrix <real> & aN ) const
89 {
90 const real xi = aXi( 0 );
91 const real eta = aXi( 1 );
92 const real zeta = aXi( 2 );
93 const real tau = 1. - xi - eta - zeta ;
94
95 real K = 64.8 ;
96
97 real eta2 = eta*eta ;
98 real zeta2 = zeta*zeta ;
99 real tau2 = tau*tau ;
100
101 real a = 2.*zeta2- eta2 -tau2 ;
102 real b = 2.*eta2 - zeta2 -tau2 ;
103 real c = 2.*tau2 - eta2 - zeta2 ;
104 real d = eta*zeta*tau ;
105
106 aN.set_size( 1, 3 );
107 aN( 0, 0 ) = K * d * a ;
108 aN( 0, 1 ) = K * d * b ;
109 aN( 0, 2 ) = K * d * c ;
110 }
111
112
113//------------------------------------------------------------------------------
114
115 template<>
116 void
119 const Vector <real> & aXi,
120 Matrix <real> & adNdXi ) const
121 {
122 const real xi = aXi( 0 );
123 const real eta = aXi( 1 );
124 const real zeta = aXi( 2 );
125 const real tau = 1. - xi - eta - zeta ;
126
127 real K = 64.8 ;
128
129 real eta2 = eta*eta ;
130 real zeta2 = zeta*zeta ;
131 real tau2 = tau*tau ;
132
133 real a = 2.*zeta2- eta2 -tau2 ;
134 real b = 2.*eta2 - zeta2 -tau2 ;
135 real c = 2.*tau2 - eta2 - zeta2 ;
136 real d = eta*zeta*tau ;
137
138 real a_xi = 2.* tau ;
139 real a_eta = 2.*( tau - eta );
140 real a_zeta = 2.* ( tau + 2 * zeta );
141
142 real b_xi = 2. * tau ;
143 real b_eta = 2. * ( tau + 2. * eta );
144 real b_zeta = 2.* ( tau - zeta ) ;
145
146 real c_xi = -4. * tau ;
147 real c_eta = -2. * ( eta + 2.*tau ) ;
148 real c_zeta = -2. * ( zeta + 2.*tau );
149
150 real d_xi = - eta * zeta ;
151 real d_eta = zeta*( tau - eta );
152 real d_zeta = eta*( tau - zeta );
153
154 adNdXi.set_size( 3, 3 );
155
156 adNdXi( 0, 0 ) = K * ( d_xi * a + d * a_xi ); // dN1/dxi
157 adNdXi( 0, 1 ) = K * ( d_eta * a + d * a_eta ); // dN1/deta
158 adNdXi( 0, 2 ) = K * ( d_zeta * a + d * a_zeta ); // dN1/dzeta
159
160 adNdXi( 1, 0 ) = K * ( d_xi * b + d * b_xi ); // dN2/dxi
161 adNdXi( 1, 1 ) = K * ( d_eta * b + d * b_eta ); // dN2/deta
162 adNdXi( 1, 2 ) = K * ( d_zeta * b + d * b_zeta ); // dN2/dzeta
163
164 adNdXi( 2, 0 ) = K * ( d_xi * c + d * c_xi ); // dN3/dxi
165 adNdXi( 2, 1 ) = K * ( d_eta * c + d * c_eta ); // dN3/deta
166 adNdXi( 2, 2 ) = K * ( d_zeta * c + d * c_zeta ); // dN3/dzeta
167
168 }
169
170//------------------------------------------------------------------------------
171
172
173 template<>
174 void
177 const Vector <real> & aXi,
178 Matrix <real> & ad2NdXi2 ) const
179 {
180 const real xi = aXi( 0 );
181 const real eta = aXi( 1 );
182 const real zeta = aXi( 2 );
183 const real tau = 1. - xi - eta - zeta ;
184
185 real K = 64.8 ;
186
187 real eta2 = eta*eta ;
188 real zeta2 = zeta*zeta ;
189 real tau2 = tau*tau ;
190
191 real a = 2.*zeta2- eta2 -tau2 ;
192 real b = 2.*eta2 - zeta2 -tau2 ;
193 real c = 2.*tau2 - eta2 - zeta2 ;
194 real d = eta*zeta*tau ;
195
196 real a_xi = 2.* tau ;
197 real a_eta = 2.*( tau - eta );
198 real a_zeta = 2.* ( tau + 2 * zeta );
199
200 real b_xi = 2. * tau ;
201 real b_eta = 2. * ( tau + 2. * eta );
202 real b_zeta = 2.* ( tau - zeta ) ;
203
204 real c_xi = -4. * tau ;
205 real c_eta = -2. * ( eta + 2.*tau ) ;
206 real c_zeta = -2. * ( zeta + 2.*tau );
207
208 real d_xi = - eta * zeta ;
209 real d_eta = zeta*( tau - eta );
210 real d_zeta = eta*( tau - zeta );
211
212 ad2NdXi2.set_size( 6,3 );
213
214 // d²N1/dxi²
215 ad2NdXi2( 0, 0 ) = 2.*K*(a_xi*d_xi-d) ;
216
217 // d²N1/deta²
218 ad2NdXi2( 1, 0 ) = 2.*K*(a_eta*d_eta - a*zeta-2.*d ) ;
219
220 // d²N1/dzeta²
221 ad2NdXi2( 2, 0 ) = 2.*K*(d + a_zeta*d_zeta - a*eta) ;
222
223
224 // d²N1/(deta*dzeta)
225 ad2NdXi2( 3, 0 ) = K*(a_eta*d_zeta + a_zeta*d_eta - a*(eta+zeta-tau)-2.*d);
226
227 // d²N1/(dxi*dzeta)
228 ad2NdXi2( 4, 0 ) = K*(a_xi*d_zeta + a_zeta*d_xi - a*eta-2.*d);
229
230 // d²N1/(dxi*deta)
231 ad2NdXi2( 5, 0 ) = K*(a_eta*d_xi + a_xi*d_eta - a*zeta-2.*d);
232
233 // d²N2/dxi²
234 ad2NdXi2( 0, 1 ) = 2.*K*(b_xi*d_xi - d );
235
236 // d²N2/deta²
237 ad2NdXi2( 1, 1 ) = 2.*K*(b_eta*d_eta - b*zeta + d);
238
239 // d²N2/dzeta²
240 ad2NdXi2( 2, 1 ) = 2.*K*(b_zeta*d_zeta - b*eta-2.*d);
241
242 // d²N2/(deta*dzeta)
243 ad2NdXi2( 3, 1 ) = K*(b_eta*d_zeta + b_zeta*d_eta - b*(eta+zeta-tau)-2.*d);
244
245 // d²N2/(dxi*dzeta)
246 ad2NdXi2( 4, 1 ) = K*(b_xi*d_zeta + b_zeta*d_xi - b*eta-2.*d);
247
248 // d²N2/(dxi*deta)
249 ad2NdXi2( 5, 1 ) = K*(b_eta*d_xi + b_xi*d_eta - b*zeta-2.*d);
250
251 // d²N3/dxi²
252 ad2NdXi2( 0, 2 ) = 2.*K*(2.*d + c_xi*d_xi);
253
254 // d²N3/deta²
255 ad2NdXi2( 1, 2 ) = 2.*K*(d + c_eta*d_eta - c*zeta);
256
257 // d²N3/dzeta²
258 ad2NdXi2( 2, 2 ) = 2.*K*(d + c_zeta*d_zeta - c*eta);
259
260 // d²N3/(deta*dzeta)
261 ad2NdXi2( 3, 2 ) = K*(4.*d + c_eta*d_zeta + c_zeta*d_eta - c*(eta+zeta-tau));
262
263 // d²N3/(dxi*dzeta)
264 ad2NdXi2( 4, 2 ) = K*(4.*d + c_xi*d_zeta + c_zeta*d_xi - c*eta);
265
266 // d²N3/(dxi*deta)
267 ad2NdXi2( 5, 2 ) = K*(4.*d + c_eta*d_xi + c_xi*d_eta - c*zeta);
268
269 }
270
271//------------------------------------------------------------------------------
272 }
273}
274#endif //BELFEM_CL_IFG_TET10B_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
@ K
Stiffness matrix.
Definition cl_TimestepMatrices.hpp:35
USER GUIDES:
Definition cl_Capacitor.cpp:16
@ BubbleFace1
Definition Mesh_Enums.hpp:107
ElementType element_type(const std::string &aStr)
Definition Mesh_Enums.hpp:370
ElementType
Element types.
Definition Mesh_Enums.hpp:27
@ TET10
Definition Mesh_Enums.hpp:39
InterpolationOrder
Definition Mesh_Enums.hpp:85
@ CUBIC
Definition Mesh_Enums.hpp:90
double real
Definition typedefs.hpp:36
@ TET
Definition Mesh_Enums.hpp:75