BELFEM 0.9.0
Berkeley Lab Finite Element Framework
Loading...
Searching...
No Matches
cl_IF_TRI15.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_TRI15_HPP
13#define BELFEM_CL_IF_TRI15_HPP
14
15
16#include "cl_IF_InterpolationFunctionTemplate.hpp"
17
18namespace belfem
19{
20 namespace fem
21 {
22//------------------------------------------------------------------------------
23
24 template<>
25 InterpolationOrder
26 InterpolationFunctionTemplate<
27 GeometryType::TRI, InterpolationType::LAGRANGE, 2, 15 >
32
33//------------------------------------------------------------------------------
34
35 template<>
36 ElementType
37 InterpolationFunctionTemplate<
38 GeometryType::TRI, InterpolationType::LAGRANGE, 2, 15 >
39 ::element_type() const
40 {
41 return ElementType::TRI15;
42 }
43
44//------------------------------------------------------------------------------
45
46 template<>
47 void
48 InterpolationFunctionTemplate<
49 GeometryType::TRI, InterpolationType::LAGRANGE, 2, 15 >
50 ::param_coords( Matrix< real > & aXiHat ) const
51 {
52 aXiHat.set_size( 2, 15 );
53
54 aXiHat( 0, 0 ) = 1.00 ;
55 aXiHat( 0, 1 ) = 0.00 ;
56 aXiHat( 0, 2 ) = 0.00 ;
57 aXiHat( 0, 3 ) = 0.75 ;
58 aXiHat( 0, 4 ) = 0.50 ;
59 aXiHat( 0, 5 ) = 0.25 ;
60 aXiHat( 0, 6 ) = 0.00 ;
61 aXiHat( 0, 7 ) = 0.00 ;
62 aXiHat( 0, 8 ) = 0.00 ;
63 aXiHat( 0, 9 ) = 0.25 ;
64 aXiHat( 0, 10 ) = 0.50 ;
65 aXiHat( 0, 11 ) = 0.75 ;
66 aXiHat( 0, 12 ) = 0.50 ;
67 aXiHat( 0, 13 ) = 0.25 ;
68 aXiHat( 0, 14 ) = 0.25 ;
69
70 aXiHat( 1, 0 ) = 0.00 ;
71 aXiHat( 1, 1 ) = 1.00 ;
72 aXiHat( 1, 2 ) = 0.00 ;
73 aXiHat( 1, 3 ) = 0.25 ;
74 aXiHat( 1, 4 ) = 0.50 ;
75 aXiHat( 1, 5 ) = 0.75 ;
76 aXiHat( 1, 6 ) = 0.75 ;
77 aXiHat( 1, 7 ) = 0.50 ;
78 aXiHat( 1, 8 ) = 0.25 ;
79 aXiHat( 1, 9 ) = 0.00 ;
80 aXiHat( 1, 10 ) = 0.00 ;
81 aXiHat( 1, 11 ) = 0.00 ;
82 aXiHat( 1, 12 ) = 0.25 ;
83 aXiHat( 1, 13 ) = 0.50 ;
84 aXiHat( 1, 14 ) = 0.25 ;
85
86 }
87
88//------------------------------------------------------------------------------
89
90 template<>
91 void
92 InterpolationFunctionTemplate<
93 GeometryType::TRI, InterpolationType::LAGRANGE, 2, 15 >::N(
94 const Vector< real > & aXi,
95 Matrix< real > & aN ) const
96 {
97 const real xi = aXi( 0 );
98 const real eta = aXi( 1 );
99
100 const real a = 1.0/3.0 ;
101 const real b = 4.0 ;
102 const real c = 16.0 / 3.0 ;
103 const real d = 32.0 ;
104
105 aN.set_size( 1, 15 );
106
107 aN( 0, 0 ) = a*xi*(2.*xi-1.)*(4.*xi-3.)*(4.*xi-1.);
108 aN( 0, 1 ) = a*eta*(2.*eta-1.)*(4.*eta-3.)*(4.*eta-1.);
109 aN( 0, 2 ) = a*(xi+eta-1.)*(2.*(xi+eta)-1.)*(4.*(xi+eta)-3.)*(4.*(xi+eta)-1.);
110 aN( 0, 3 ) = c*eta*xi*(8.*xi*xi-6.*xi+1.);
111 aN( 0, 4 ) = b*eta*xi*(4.*eta-1.)*(4.*xi-1.);
112 aN( 0, 5 ) = c*eta*xi*(8.*eta*eta-6.*eta+1.);
113 aN( 0, 6 ) = -c*eta*(xi+eta-1.)*(8.*eta*eta-6.*eta+1.);
114 aN( 0, 7 ) = b*eta*(4.*eta-1.)*(xi+eta-1.)*(4.*(xi+eta)-3.);
115 aN( 0, 8 ) = -c*eta*(xi+eta-1.)*(2.*(xi+eta)-1.)*(4.*(xi+eta)-3.);
116 aN( 0, 9 ) = -c*xi*(xi+eta-1.)*(2.*(xi+eta)-1.)*(4.*(xi+eta)-3.);
117 aN( 0, 10 ) = b*xi*(xi+eta-1.)*(4.*xi-1.)*(4.*(xi+eta)-3.);
118 aN( 0, 11 ) = -c*xi*(xi+eta-1.)*(8.*xi*xi-6.*xi+1.);
119 aN( 0, 12 ) = -d*eta*xi*(4.*xi-1.)*(xi+eta-1.);
120 aN( 0, 13 ) = -d*eta*xi*(4.*eta-1.)*(xi+eta-1.);
121 aN( 0, 14 ) = d*eta*xi*(xi+eta-1.)*(4.*(xi+eta)-3.);
122 }
123
124//------------------------------------------------------------------------------
125
126 template<>
127 void
128 InterpolationFunctionTemplate<
129 GeometryType::TRI, InterpolationType::LAGRANGE, 2, 15 >::dNdXi(
130 const Vector< real > & aXi,
131 Matrix< real > & adNdXi ) const
132 {
133 adNdXi.set_size( 2, 15 );
134
135 const real xi = aXi( 0 );
136 const real eta = aXi( 1 );
137
138 const real a = 1.0/3.0 ;
139 const real b = 4.0 ;
140 const real c = 16.0 / 3.0 ;
141 const real d = 32.0 ;
142
143 adNdXi( 0, 0 ) = a*(8.*xi-3.)*(1.+4.*xi*(4.*xi-3.)) ;
144 adNdXi( 1, 0 ) = 0.0 ;
145
146 adNdXi( 0, 1 ) = 0.0 ;
147 adNdXi( 1, 1 ) = a*(8.*eta-3.)*(1+4.*eta*(4.*eta-3.)) ;
148
149 adNdXi( 0, 2 ) = a*(8.*(xi+eta)-5.)*(5.+16.*eta*eta+4.*xi*(4.*xi-5.)+4.*eta*(8.*xi-5.)) ;
150 adNdXi( 1, 2 ) = a*(8.*(xi+eta)-5.)*(5.+4*xi*(4.*xi-5.)+4*eta*((8.*xi-5.)+4.*eta)) ;
151
152 adNdXi( 0, 3 ) = c*eta*(1.+12.*xi*(2.*xi-1.)) ;
153 adNdXi( 1, 3 ) = c*xi*(8.*xi*xi-6.*xi+1.) ;
154
155 adNdXi( 0, 4 ) = b*eta*(4.*eta-1.)*(8.*xi-1.) ;
156 adNdXi( 1, 4 ) = b*xi*(8.*eta-1.)*(4.*xi-1.) ;
157
158 adNdXi( 0, 5 ) = c*eta*(eta*(8.*eta-6.)+1.) ;
159 adNdXi( 1, 5 ) = c*xi*(24.*eta*eta-12*eta+1.) ;
160
161 adNdXi( 0, 6 ) = -c*eta*(eta*(8.*eta-6.)+1.) ;
162 adNdXi( 1, 6 ) = -(c*(xi+2.*eta*(7.-6.*xi+eta*(16*eta+12.*xi-21.))-1.)) ;
163
164 adNdXi( 0, 7 ) = b*eta*(4.*eta-1.)*(8.*(eta+xi)-7.) ;
165 adNdXi( 1, 7 ) = b*(2.*eta+xi-1.)*(3.-4.*xi+32*eta*(xi+eta-1.)) ;
166
167 adNdXi( 0, 8 ) = -(c*eta*(13.+24.*eta*eta+12.*xi*(2.*xi-3.)+12.*eta*(4.*xi-3.))) ;
168 adNdXi( 1, 8 ) = -(c*(-3.+2*eta*(13.+eta*(16*eta-27.))+xi*(13.+72.*(eta-1.)*eta+xi*(6.*(8.*eta-3.)+8.*xi)))) ;
169
170 adNdXi( 0, 9 ) = -(c*(eta*(eta*(8.*eta+6.*(8.*xi-3.))+(13.+72.*(xi-1.)*xi))+2.*xi*(13.+xi*(16.*xi-27.))-3.)) ;
171 adNdXi( 1, 9 ) = -(c*xi*(13+24.*eta*eta+12.*xi*(2*xi-3.)+12.*eta*(4.*xi-3.))) ;
172
173 adNdXi( 0, 10 ) = b*(2.*xi+eta-1.)*(3.+32.*(xi-1.)*xi+4.*eta*(8.*xi-1.)) ;
174 adNdXi( 1, 10 ) = b*xi*(4.*xi-1.)*(8.*eta+8.*xi-7.) ;
175
176 adNdXi( 0, 11 ) = -(c*(eta+12.*eta*xi*(2.*xi-1.)+2.*xi*(7.+xi*(16.*xi-21.))-1.)) ;
177 adNdXi( 1, 11 ) = -c*xi*(8.*xi*xi-6.*xi+1.) ;
178
179 adNdXi( 0, 12 ) = d*eta*(eta-8.*eta*xi+2.*(5.-6.*xi)*xi-1.) ;
180 adNdXi( 1, 12 ) = -d*xi*(4.*xi-1.)*(2*eta+xi-1.) ;
181
182 adNdXi( 0, 13 ) = -(d*eta*(4.*eta-1.)*(2.*xi+eta-1.)) ;
183 adNdXi( 1, 13 ) = d*xi*(xi-2*eta*(6*eta+4.*xi-5.)-1.) ;
184
185 adNdXi( 0, 14 ) = d*eta*(3.+4.*eta*eta+2.*xi*(6.*xi-7.)+eta*(16.*xi-7.)) ;
186 adNdXi( 1, 14 ) = d*xi*(3.+2*eta*(6*eta-7.)+xi*(4.*xi-7.+16.*eta)) ;
187 }
188
189//------------------------------------------------------------------------------
190
191 template<>
192 void
193 InterpolationFunctionTemplate<
194 GeometryType::TRI, InterpolationType::LAGRANGE, 2, 15 >::d2NdXi2(
195 const Vector< real > & aXi,
196 Matrix< real > & ad2NdXi2 ) const
197 {
198 ad2NdXi2.set_size( 3, 15, 0.0 );
199
200 const real xi = aXi( 0 );
201 const real eta = aXi( 1 );
202
203 const real a = 1.0/3.0 ;
204 const real b = 4.0 ;
205 const real c = 16.0 / 3.0 ;
206 const real d = 32.0 ;
207
208 ad2NdXi2( 0, 0 ) = 4.*a*(11.+24.*xi*(4.*xi-3.)) ;
209 ad2NdXi2( 1, 0 ) = 0.0 ;
210 ad2NdXi2( 2, 0 ) = 0.0 ;
211
212 ad2NdXi2( 0, 1 ) = 0.0 ;
213 ad2NdXi2( 1, 1 ) = 4.*a*(11.+24.*eta*(4.*eta-3.)) ;
214 ad2NdXi2( 2, 1 ) = 0.0 ;
215
216 ad2NdXi2( 0, 2 ) = 4.*a*(35.+24.*xi*(4.*xi-5.)+eta*(96.*eta+24*(8.*xi-5.))) ;
217 ad2NdXi2( 1, 2 ) = 4.*a*(35.+24.*xi*(4.*xi-5.)+eta*(96.*eta+24*(8.*xi-5.))) ;
218 ad2NdXi2( 2, 2 ) = 4.*a*(35.+24*xi*(4.*xi-5.)+24.*eta*(4.*eta+8.*xi-5.)) ;
219
220 ad2NdXi2( 0, 3 ) = c*eta*(48.*xi-12.) ;
221 ad2NdXi2( 1, 3 ) = 0.0 ;
222 ad2NdXi2( 2, 3 ) = c*(12.*xi*(2.*xi-1)+1.) ;
223
224 ad2NdXi2( 0, 4 ) = 8.*b*eta*(4.*eta-1.) ;
225 ad2NdXi2( 1, 4 ) = 8.*b*xi*(4.*xi-1.) ;
226 ad2NdXi2( 2, 4 ) = b*(8.*eta-1.)*(8.*xi-1.) ;
227
228 ad2NdXi2( 0, 5 ) = 0.0;
229 ad2NdXi2( 1, 5 ) = c*xi*(48.*eta-12.) ;
230 ad2NdXi2( 2, 5 ) = c+12.*c*eta*(2.*eta-1.) ;
231
232 ad2NdXi2( 0, 6 ) = 0.0;
233 ad2NdXi2( 1, 6 ) = 2.*c*(6*xi-6.*eta*(8.*eta+4.*xi-7.)-7) ;
234 ad2NdXi2( 2, 6 ) = c*(12.*(1.-2.*eta)*eta-1.) ;
235
236 ad2NdXi2( 0, 7 ) = 8.*b*eta*(4.*eta-1.) ;
237 ad2NdXi2( 1, 7 ) = 2.*b*(19.+96.*eta*(xi+eta-1.)+4.*xi*(4.*xi-9.)) ;
238 ad2NdXi2( 2, 7 ) = b*(7.-8.*xi+8.*eta*(12.*eta+8.*xi-9.)) ;
239
240 ad2NdXi2( 0, 8 ) = c*eta*(36.-48.*(xi+eta)) ;
241 ad2NdXi2( 1, 8 ) = 2.*c*(36.*xi-13.-6.*eta*(8.*eta-9.)-24.*xi*(3.*eta+xi)) ;
242 ad2NdXi2( 2, 8 ) = c*(36.*xi-24.*(3.*eta*(eta-1.)+xi*(4.*eta+xi))-13.) ;
243
244 ad2NdXi2( 0, 9 ) = 2.*c*(6.*xi*(9.-8.*xi)-eta*(24.*eta+36.*(2.*xi-1.))-13.) ;
245 ad2NdXi2( 1, 9 ) = c*xi*(36-48*(xi+eta)) ;
246 ad2NdXi2( 2, 9 ) = c*(72.*(1.-xi)*xi-12.*eta*(2.*eta+8.*xi-3.)-13.) ;
247
248 ad2NdXi2( 0, 10 ) = 2.*b*(19.+96.*(xi-1.)*xi+eta*(16.*eta+12.*(8.*xi-3.))) ;
249 ad2NdXi2( 1, 10 ) = 8.*b*xi*(4.*xi-1.) ;
250 ad2NdXi2( 2, 10 ) = b*(7.+24.*xi*(4.*xi-3.)+8.*eta*(8.*xi-1.)) ;
251
252 ad2NdXi2( 0, 11 ) = 2.*c*(6.*eta*(1.0-4.*xi)-6.*xi*(8.*xi-7.)-7.) ;
253 ad2NdXi2( 1, 11 ) = 0.0 ;
254 ad2NdXi2( 2, 11 ) = c*(12.*xi*(1.-2.*xi)-1.) ;
255
256 ad2NdXi2( 0, 12 ) = d*eta*(10-24.*xi-8.*eta) ;
257 ad2NdXi2( 1, 12 ) = 2.*d*xi*(1.0-4.*xi) ;
258 ad2NdXi2( 2, 12 ) = d*(eta*(2.-16.*xi)+2.*(5. - 6.*xi)*xi -1.) ;
259
260 ad2NdXi2( 0, 13 ) = 2.*d*eta*(1.-4.*eta) ;
261 ad2NdXi2( 1, 13 ) = d*xi*(10.-24.*eta-8.*xi) ;
262 ad2NdXi2( 2, 13 ) = d*(2.*xi - 2.*eta*(6.*eta+8.*xi-5.)-1.) ;
263
264 ad2NdXi2( 0, 14 ) = d*eta*(16.*eta+24.*xi-14.) ;
265 ad2NdXi2( 1, 14 ) = d*xi*(24.*eta+16.*xi-14.) ;
266 ad2NdXi2( 2, 14 ) = d*(3.+2*xi*(6.*xi-7.) + 2*eta*(6.*eta+16*xi-7.)) ;
267 }
268
269//------------------------------------------------------------------------------
270 }
271}
272
273#endif //BELFEM_CL_IF_TRI15_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::TRI15
@ TRI15
Definition Mesh_Enums.hpp:49
belfem::InterpolationOrder
InterpolationOrder
Definition Mesh_Enums.hpp:85
belfem::InterpolationOrder::QUARTIC
@ QUARTIC
Definition Mesh_Enums.hpp:91
belfem::real
double real
Definition typedefs.hpp:36
belfem::GeometryType::TRI
@ TRI
Definition Mesh_Enums.hpp:73