svMultiPhysics
Loading...
Searching...
No Matches
nn_elem_gnn.h
1// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the University of California, and others.
2// SPDX-License-Identifier: BSD-3-Clause
3
4/// @brief Define a map type used to set element shape function data.
5///
6/// Reproduces the Fortran 'GETGNN' subroutine.
7//
8using GetElementShapeMapType = std::map<ElementType, std::function<void(const int, const int, const int,
9 Array<double>&, Array<double>&, Array3<double>&)>>;
10
11GetElementShapeMapType get_element_shape_data = {
12
13 {ElementType::HEX8, [](const int insd, const int eNoN, const int g, Array<double>& xi, Array<double>& N,
14 Array3<double>& Nx) -> void
15 {
16 double lx = 1.0 - xi(0,g);
17 double ly = 1.0 - xi(1,g);
18 double lz = 1.0 - xi(2,g);
19 double ux = 1.0 + xi(0,g);
20 double uy = 1.0 + xi(1,g);
21 double uz = 1.0 + xi(2,g);
22
23 N(0,g) = lx*ly*lz/8.0;
24 N(1,g) = ux*ly*lz/8.0;
25 N(2,g) = ux*uy*lz/8.0;
26 N(3,g) = lx*uy*lz/8.0;
27 N(4,g) = lx*ly*uz/8.0;
28 N(5,g) = ux*ly*uz/8.0;
29 N(6,g) = ux*uy*uz/8.0;
30 N(7,g) = lx*uy*uz/8.0;
31
32 Nx(0,0,g) = -ly*lz/8.0;
33 Nx(1,0,g) = -lx*lz/8.0;
34 Nx(2,0,g) = -lx*ly/8.0;
35
36 Nx(0,1,g) = ly*lz/8.0;
37 Nx(1,1,g) = -ux*lz/8.0;
38 Nx(2,1,g) = -ux*ly/8.0;
39
40 Nx(0,2,g) = uy*lz/8.0;
41 Nx(1,2,g) = ux*lz/8.0;
42 Nx(2,2,g) = -ux*uy/8.0;
43
44 Nx(0,3,g) = -uy*lz/8.0;
45 Nx(1,3,g) = lx*lz/8.0;
46 Nx(2,3,g) = -lx*uy/8.0;
47
48 Nx(0,4,g) = -ly*uz/8.0;
49 Nx(1,4,g) = -lx*uz/8.0;
50 Nx(2,4,g) = lx*ly/8.0;
51
52 Nx(0,5,g) = ly*uz/8.0;
53 Nx(1,5,g) = -ux*uz/8.0;
54 Nx(2,5,g) = ux*ly/8.0;
55
56 Nx(0,6,g) = uy*uz/8.0;
57 Nx(1,6,g) = ux*uz/8.0;
58 Nx(2,6,g) = ux*uy/8.0;
59
60 Nx(0,7,g) = -uy*uz/8.0;
61 Nx(1,7,g) = lx*uz/8.0;
62 Nx(2,7,g) = lx*uy/8.0;
63 }
64 },
65
66 {ElementType::HEX20, [](const int insd, const int eNoN, const int g, Array<double>& xi, Array<double>& N,
67 Array3<double>& Nx) -> void
68 {
69 double lx = 1.0 - xi(0,g);
70 double ly = 1.0 - xi(1,g);
71 double lz = 1.0 - xi(2,g);
72 double ux = 1.0 + xi(0,g);
73 double uy = 1.0 + xi(1,g);
74 double uz = 1.0 + xi(2,g);
75
76 double mx = lx*ux;
77 double my = ly*uy;
78 double mz = lz*uz;
79
80 N(0, g) = lx*ly*lz*(lx+ly+lz-5.0)/8.0;
81 N(1, g) = ux*ly*lz*(ux+ly+lz-5.0)/8.0;
82 N(2, g) = ux*uy*lz*(ux+uy+lz-5.0)/8.0;
83 N(3, g) = lx*uy*lz*(lx+uy+lz-5.0)/8.0;
84 N(4, g) = lx*ly*uz*(lx+ly+uz-5.0)/8.0;
85 N(5, g) = ux*ly*uz*(ux+ly+uz-5.0)/8.0;
86 N(6, g) = ux*uy*uz*(ux+uy+uz-5.0)/8.0;
87 N(7, g) = lx*uy*uz*(lx+uy+uz-5.0)/8.0;
88 N(8, g) = mx*ly*lz/4.0;
89 N(9, g) = ux*my*lz/4.0;
90 N(10, g) = mx*uy*lz/4.0;
91 N(11, g) = lx*my*lz/4.0;
92 N(12, g) = mx*ly*uz/4.0;
93 N(13, g) = ux*my*uz/4.0;
94 N(14, g) = mx*uy*uz/4.0;
95 N(15, g) = lx*my*uz/4.0;
96 N(16, g) = lx*ly*mz/4.0;
97 N(17, g) = ux*ly*mz/4.0;
98 N(18, g) = ux*uy*mz/4.0;
99 N(19, g) = lx*uy*mz/4.0;
100
101 // N(1) = lx*ly*lz*(lx+ly+lz-5.0)/8.0;
102 int n = 0;
103 Nx(0,n,g) = -ly*lz*(lx+ly+lz-5.0+lx)/8.0;
104 Nx(1,n,g) = -lx*lz*(lx+ly+lz-5.0+ly)/8.0;
105 Nx(2,n,g) = -lx*ly*(lx+ly+lz-5.0+lz)/8.0;
106
107//c N(n,g) = ux*ly*lz*(ux+ly+lz-5.0)/8.0;
108 n += 1;
109 Nx(0,n,g) = ly*lz*(ux+ly+lz-5.0+ux)/8.0;
110 Nx(1,n,g) = -ux*lz*(ux+ly+lz-5.0+ly)/8.0;
111 Nx(2,n,g) = -ux*ly*(ux+ly+lz-5.0+lz)/8.0;
112
113//c N(n,g) = ux*uy*lz*(ux+uy+lz-5.0)/8.0
114 n += 1;
115 Nx(0,n,g) = uy*lz*(ux+uy+lz-5.0+ux)/8.0;
116 Nx(1,n,g) = ux*lz*(ux+uy+lz-5.0+uy)/8.0;
117 Nx(2,n,g) = -ux*uy*(ux+uy+lz-5.0+lz)/8.0;
118
119//c N(n,g) = lx*uy*lz*(lx+uy+lz-5.0)/8.0
120 n += 1;
121 Nx(0,n,g) = -uy*lz*(lx+uy+lz-5.0+lx)/8.0;
122 Nx(1,n,g) = lx*lz*(lx+uy+lz-5.0+uy)/8.0;
123 Nx(2,n,g) = -lx*uy*(lx+uy+lz-5.0+lz)/8.0;
124
125//c N(n,g) = lx*ly*uz*(lx+ly+uz-5.0)/8.0
126 n += 1;
127 Nx(0,n,g) = -ly*uz*(lx+ly+uz-5.0+lx)/8.0;
128 Nx(1,n,g) = -lx*uz*(lx+ly+uz-5.0+ly)/8.0;
129 Nx(2,n,g) = lx*ly*(lx+ly+uz-5.0+uz)/8.0;
130
131//c N(n,g) = ux*ly*uz*(ux+ly+uz-5.0)/8.0
132 n += 1;
133 Nx(0,n,g) = ly*uz*(ux+ly+uz-5.0+ux)/8.0;
134 Nx(1,n,g) = -ux*uz*(ux+ly+uz-5.0+ly)/8.0;
135 Nx(2,n,g) = ux*ly*(ux+ly+uz-5.0+uz)/8.0;
136
137//c N(n,g) = ux*uy*uz*(ux+uy+uz-5.0)/8.0
138 n += 1;
139 Nx(0,n,g) = uy*uz*(ux+uy+uz-5.0+ux)/8.0;
140 Nx(1,n,g) = ux*uz*(ux+uy+uz-5.0+uy)/8.0;
141 Nx(2,n,g) = ux*uy*(ux+uy+uz-5.0+uz)/8.0;
142
143//c N(n,g) = lx*uy*uz*(lx+uy+uz-5.0)/8.0
144 n += 1;
145 Nx(0,n,g) = -uy*uz*(lx+uy+uz-5.0+lx)/8.0;
146 Nx(1,n,g) = lx*uz*(lx+uy+uz-5.0+uy)/8.0;
147 Nx(2,n,g) = lx*uy*(lx+uy+uz-5.0+uz)/8.0;
148
149//c N(n,g) = mx*ly*lz/4.0
150 n += 1;
151 Nx(0,n,g) = (lx - ux)*ly*lz/4.0;
152 Nx(1,n,g) = -mx*lz/4.0;
153 Nx(2,n,g) = -mx*ly/4.0;
154
155//c N(0n,g) = ux*my*lz/4.0
156 n += 1;
157 Nx(0,n,g) = my*lz/4.0;
158 Nx(1,n,g) = (ly - uy)*ux*lz/4.0;
159 Nx(2,n,g) = -ux*my/4.0;
160
161//c N(0n,g) = mx*uy*lz/4.0
162 n += 1;
163 Nx(0,n,g) = (lx - ux)*uy*lz/4.0;
164 Nx(1,n,g) = mx*lz/4.0;
165 Nx(2,n,g) = -mx*uy/4.0;
166
167//c N(0n,g) = lx*my*lz/4.0
168 n += 1;
169 Nx(0,n,g) = -my*lz/4.0;
170 Nx(1,n,g) = (ly - uy)*lx*lz/4.0;
171 Nx(2,n,g) = -lx*my/4.0;
172
173//c N(0n,g) = mx*ly*uz/4.0
174 n += 1;
175 Nx(0,n,g) = (lx - ux)*ly*uz/4.0;
176 Nx(1,n,g) = -mx*uz/4.0;
177 Nx(2,n,g) = mx*ly/4.0;
178
179//c N(0n,g) = ux*my*uz/4.0
180 n += 1;
181 Nx(0,n,g) = my*uz/4.0;
182 Nx(1,n,g) = (ly - uy)*ux*uz/4.0;
183 Nx(2,n,g) = ux*my/4.0;
184
185//c N(0n,g) = mx*uy*uz/4.0
186 n += 1;
187 Nx(0,n,g) = (lx - ux)*uy*uz/4.0;
188 Nx(1,n,g) = mx*uz/4.0;
189 Nx(2,n,g) = mx*uy/4.0;
190
191//c N(0n,g) = lx*my*uz/4.0
192 n += 1;
193 Nx(0,n,g) = -my*uz/4.0;
194 Nx(1,n,g) = (ly - uy)*lx*uz/4.0;
195 Nx(2,n,g) = lx*my/4.0;
196
197//c N(0n,g) = lx*ly*mz/4.0
198 n += 1;
199 Nx(0,n,g) = -ly*mz/4.0;
200 Nx(1,n,g) = -lx*mz/4.0;
201 Nx(2,n,g) = (lz - uz)*lx*ly/4.0;
202
203//c N(0n,g) = ux*ly*mz/4.0
204 n += 1;
205 Nx(0,n,g) = ly*mz/4.0;
206 Nx(1,n,g) = -ux*mz/4.0;
207 Nx(2,n,g) = (lz - uz)*ux*ly/4.0;
208
209//c N(0n,g) = ux*uy*mz/4.0
210 n += 1;
211 Nx(0,n,g) = uy*mz/4.0;
212 Nx(1,n,g) = ux*mz/4.0;
213 Nx(2,n,g) = (lz - uz)*ux*uy/4.0;
214
215//c N(n,g) = lx*uy*mz/4.0
216 n += 1;
217 Nx(0,n,g) = -uy*mz/4.0;
218 Nx(1,n,g) = lx*mz/4.0;
219 Nx(2,n,g) = (lz - uz)*lx*uy/4.0;
220 }
221 },
222
223 {ElementType::HEX27, [](const int insd, const int eNoN, const int g, Array<double>& xi, Array<double>& N,
224 Array3<double>& Nx) -> void
225 {
226 double lx = 1.0 - xi(0,g);
227 double ly = 1.0 - xi(1,g);
228 double lz = 1.0 - xi(2,g);
229 double ux = 1.0 + xi(0,g);
230 double uy = 1.0 + xi(1,g);
231 double uz = 1.0 + xi(2,g);
232
233 double mx = xi(0,g);
234 double my = xi(1,g);
235 double mz = xi(2,g);
236
237 N(0,g) = -mx*lx*my*ly*mz*lz/8.0;
238 N(1,g) = mx*ux*my*ly*mz*lz/8.0;
239 N(2,g) = -mx*ux*my*uy*mz*lz/8.0;
240 N(3,g) = mx*lx*my*uy*mz*lz/8.0;
241 N(4,g) = mx*lx*my*ly*mz*uz/8.0;
242 N(5,g) = -mx*ux*my*ly*mz*uz/8.0;
243 N(6,g) = mx*ux*my*uy*mz*uz/8.0;
244 N(7,g) = -mx*lx*my*uy*mz*uz/8.0;
245 N(8,g) = lx*ux*my*ly*mz*lz/4.0;
246 N(9,g) = -mx*ux*ly*uy*mz*lz/4.0;
247 N(10,g) = -lx*ux*my*uy*mz*lz/4.0;
248 N(11,g) = mx*lx*ly*uy*mz*lz/4.0;
249 N(12,g) = -lx*ux*my*ly*mz*uz/4.0;
250 N(13,g) = mx*ux*ly*uy*mz*uz/4.0;
251 N(14,g) = lx*ux*my*uy*mz*uz/4.0;
252 N(15,g) = -mx*lx*ly*uy*mz*uz/4.0;
253 N(16,g) = mx*lx*my*ly*lz*uz/4.0;
254 N(17,g) = -mx*ux*my*ly*lz*uz/4.0;
255 N(18,g) = mx*ux*my*uy*lz*uz/4.0;
256 N(19,g) = -mx*lx*my*uy*lz*uz/4.0;
257
258 N(20,g) = -mx*lx*ly*uy*lz*uz/2.0;
259 N(21,g) = mx*ux*ly*uy*lz*uz/2.0;
260 N(22,g) = -lx*ux*my*ly*lz*uz/2.0;
261 N(23,g) = lx*ux*my*uy*lz*uz/2.0;
262 N(24,g) = -lx*ux*ly*uy*mz*lz/2.0;
263 N(25,g) = lx*ux*ly*uy*mz*uz/2.0;
264
265 N(26,g) = lx*ux*ly*uy*lz*uz;
266
267 // N(0) = -mx*lx*my*ly*mz*lz/8.0
268 int n = 0;
269 Nx(0,n,g) = -(lx - mx)*my*ly*mz*lz/8.0;
270 Nx(1,n,g) = -(ly - my)*mx*lx*mz*lz/8.0;
271 Nx(2,n,g) = -(lz - mz)*mx*lx*my*ly/8.0;
272
273 // N(n,g) = mx*ux*my*ly*mz*lz/8.0
274 n += 1;
275 Nx(0,n,g) = (mx + ux)*my*ly*mz*lz/8.0;
276 Nx(1,n,g) = (ly - my)*mx*ux*mz*lz/8.0;
277 Nx(2,n,g) = (lz - mz)*mx*ux*my*ly/8.0;
278
279 // N(n,g) = -mx*ux*my*uy*mz*lz/8.0
280 n += 1;
281 Nx(0,n,g) = -(mx + ux)*my*uy*mz*lz/8.0;
282 Nx(1,n,g) = -(my + uy)*mx*ux*mz*lz/8.0;
283 Nx(2,n,g) = -(lz - mz)*mx*ux*my*uy/8.0;
284
285 // N(n,g) = mx*lx*my*uy*mz*lz/8.0
286 n += 1;
287 Nx(0,n,g) = (lx - mx)*my*uy*mz*lz/8.0;
288 Nx(1,n,g) = (my + uy)*mx*lx*mz*lz/8.0;
289 Nx(2,n,g) = (lz - mz)*mx*lx*my*uy/8.0;
290
291 // N(n,g) = mx*lx*my*ly*mz*uz/8.0
292 n += 1;
293 Nx(0,n,g) = (lx - mx)*my*ly*mz*uz/8.0;
294 Nx(1,n,g) = (ly - my)*mx*lx*mz*uz/8.0;
295 Nx(2,n,g) = (mz + uz)*mx*lx*my*ly/8.0;
296
297 // N(n,g) = -mx*ux*my*ly*mz*uz/8.0
298 n += 1;
299 Nx(0,n,g) = -(mx + ux)*my*ly*mz*uz/8.0;
300 Nx(1,n,g) = -(ly - my)*mx*ux*mz*uz/8.0;
301 Nx(2,n,g) = -(mz + uz)*mx*ux*my*ly/8.0;
302
303 // N(n,g) = mx*ux*my*uy*mz*uz/8.0
304 n += 1;
305 Nx(0,n,g) = (mx + ux)*my*uy*mz*uz/8.0;
306 Nx(1,n,g) = (my + uy)*mx*ux*mz*uz/8.0;
307 Nx(2,n,g) = (mz + uz)*mx*ux*my*uy/8.0;
308
309 // N(n,g) = -mx*lx*my*uy*mz*uz/8.0
310 n += 1;
311 Nx(0,n,g) = -(lx - mx)*my*uy*mz*uz/8.0;
312 Nx(1,n,g) = -(my + uy)*mx*lx*mz*uz/8.0;
313 Nx(2,n,g) = -(mz + uz)*mx*lx*my*uy/8.0;
314
315 // N(n,g) = lx*ux*my*ly*mz*lz/4.0
316 n += 1;
317 Nx(0,n,g) = (lx - ux)*my*ly*mz*lz/4.0;
318 Nx(1,n,g) = (ly - my)*lx*ux*mz*lz/4.0;
319 Nx(2,n,g) = (lz - mz)*lx*ux*my*ly/4.0;
320
321 // N(n,g) = -mx*ux*ly*uy*mz*lz/4.0
322 n += 1;
323 Nx(0,n,g) = -(mx + ux)*ly*uy*mz*lz/4.0;
324 Nx(1,n,g) = -(ly - uy)*mx*ux*mz*lz/4.0;
325 Nx(2,n,g) = -(lz - mz)*mx*ux*ly*uy/4.0;
326
327 // N(n,g) = -lx*ux*my*uy*mz*lz/4.0
328 n += 1;
329 Nx(0,n,g) = -(lx - ux)*my*uy*mz*lz/4.0;
330 Nx(1,n,g) = -(my + uy)*lx*ux*mz*lz/4.0;
331 Nx(2,n,g) = -(lz - mz)*lx*ux*my*uy/4.0;
332
333 // N(n,g) = mx*lx*ly*uy*mz*lz/4.0
334 n += 1;
335 Nx(0,n,g) = (lx - mx)*ly*uy*mz*lz/4.0;
336 Nx(1,n,g) = (ly - uy)*mx*lx*mz*lz/4.0;
337 Nx(2,n,g) = (lz - mz)*mx*lx*ly*uy/4.0;
338
339 // N(n,g) = -lx*ux*my*ly*mz*uz/4.0
340 n += 1;
341 Nx(0,n,g) = -(lx - ux)*my*ly*mz*uz/4.0;
342 Nx(1,n,g) = -(ly - my)*lx*ux*mz*uz/4.0;
343 Nx(2,n,g) = -(mz + uz)*lx*ux*my*ly/4.0;
344
345 // N(n,g) = mx*ux*ly*uy*mz*uz/4.0
346 n += 1;
347 Nx(0,n,g) = (mx + ux)*ly*uy*mz*uz/4.0;
348 Nx(1,n,g) = (ly - uy)*mx*ux*mz*uz/4.0;
349 Nx(2,n,g) = (mz + uz)*mx*ux*ly*uy/4.0;
350
351 // N(n,g) = lx*ux*my*uy*mz*uz/4.0
352 n += 1;
353 Nx(0,n,g) = (lx - ux)*my*uy*mz*uz/4.0;
354 Nx(1,n,g) = (my + uy)*lx*ux*mz*uz/4.0;
355 Nx(2,n,g) = (mz + uz)*lx*ux*my*uy/4.0;
356
357 // N(n,g) = -mx*lx*ly*uy*mz*uz/4.0
358 n += 1;
359 Nx(0,n,g) = -(lx - mx)*ly*uy*mz*uz/4.0;
360 Nx(1,n,g) = -(ly - uy)*mx*lx*mz*uz/4.0;
361 Nx(2,n,g) = -(mz + uz)*mx*lx*ly*uy/4.0;
362
363 // N(n,g) = mx*lx*my*ly*lz*uz/4.0
364 n += 1;
365 Nx(0,n,g) = (lx - mx)*my*ly*lz*uz/4.0;
366 Nx(1,n,g) = (ly - my)*mx*lx*lz*uz/4.0;
367 Nx(2,n,g) = (lz - uz)*mx*lx*my*ly/4.0;
368
369 // N(n,g) = -mx*ux*my*ly*lz*uz/4.0
370 n += 1;
371 Nx(0,n,g) = -(mx + ux)*my*ly*lz*uz/4.0;
372 Nx(1,n,g) = -(ly - my)*mx*ux*lz*uz/4.0;
373 Nx(2,n,g) = -(lz - uz)*mx*ux*my*ly/4.0;
374
375 // N(n,g) = mx*ux*my*uy*lz*uz/4.0
376 n += 1;
377 Nx(0,n,g) = (mx + ux)*my*uy*lz*uz/4.0;
378 Nx(1,n,g) = (my + uy)*mx*ux*lz*uz/4.0;
379 Nx(2,n,g) = (lz - uz)*mx*ux*my*uy/4.0;
380
381 // N(n,g) = -mx*lx*my*uy*lz*uz/4.0
382 n += 1;
383 Nx(0,n,g) = -(lx - mx)*my*uy*lz*uz/4.0;
384 Nx(1,n,g) = -(my + uy)*mx*lx*lz*uz/4.0;
385 Nx(2,n,g) = -(lz - uz)*mx*lx*my*uy/4.0;
386
387 // N(n,g) = -mx*lx*ly*uy*lz*uz/2.0
388 n += 1;
389 Nx(0,n,g) = -(lx - mx)*ly*uy*lz*uz/2.0;
390 Nx(1,n,g) = -(ly - uy)*mx*lx*lz*uz/2.0;
391 Nx(2,n,g) = -(lz - uz)*mx*lx*ly*uy/2.0;
392
393 // N(n,g) = mx*ux*ly*uy*lz*uz/2.0
394 n += 1;
395 Nx(0,n,g) = (mx + ux)*ly*uy*lz*uz/2.0;
396 Nx(1,n,g) = (ly - uy)*mx*ux*lz*uz/2.0;
397 Nx(2,n,g) = (lz - uz)*mx*ux*ly*uy/2.0;
398
399 // N(n,g) = -lx*ux*my*ly*lz*uz/2.0
400 n += 1;
401 Nx(0,n,g) = -(lx - ux)*my*ly*lz*uz/2.0;
402 Nx(1,n,g) = -(ly - my)*lx*ux*lz*uz/2.0;
403 Nx(2,n,g) = -(lz - uz)*lx*ux*my*ly/2.0;
404
405 // N(n,g) = lx*ux*my*uy*lz*uz/2.0
406 n += 1;
407 Nx(0,n,g) = (lx - ux)*my*uy*lz*uz/2.0;
408 Nx(1,n,g) = (my + uy)*lx*ux*lz*uz/2.0;
409 Nx(2,n,g) = (lz - uz)*lx*ux*my*uy/2.0;
410
411 // N(n,g) = -lx*ux*ly*uy*mz*lz/2.0
412 n += 1;
413 Nx(0,n,g) = -(lx - ux)*ly*uy*mz*lz/2.0;
414 Nx(1,n,g) = -(ly - uy)*lx*ux*mz*lz/2.0;
415 Nx(2,n,g) = -(lz - mz)*lx*ux*ly*uy/2.0;
416
417 // N(n,g) = lx*ux*ly*uy*mz*uz/2.0
418 n += 1;
419 Nx(0,n,g) = (lx - ux)*ly*uy*mz*uz/2.0;
420 Nx(1,n,g) = (ly - uy)*lx*ux*mz*uz/2.0;
421 Nx(2,n,g) = (mz + uz)*lx*ux*ly*uy/2.0;
422
423 // N(n,g) = lx*ux*ly*uy*lz*uz
424 n += 1;
425 Nx(0,n,g) = (lx - ux)*ly*uy*lz*uz;
426 Nx(1,n,g) = (ly - uy)*lx*ux*lz*uz;
427 Nx(2,n,g) = (lz - uz)*lx*ux*ly*uy;
428 }
429 },
430
431 {ElementType::LIN1, [](const int insd, const int eNoN, const int g, Array<double>& xi, Array<double>& N,
432 Array3<double>& Nx) -> void
433 {
434 N(0,g) = (1.0 - xi(0,g))*0.5;
435 N(1,g) = (1.0 + xi(0,g))*0.5;
436
437 Nx(0,0,g) = -0.5;
438 Nx(0,1,g) = 0.5;
439 }
440 },
441
442 {ElementType::LIN2, [](const int insd, const int eNoN, const int g, Array<double>& xi, Array<double>& N,
443 Array3<double>& Nx) -> void
444 {
445 N(0,g) = -xi(0,g)*(1.0 - xi(0,g))*0.50;
446 N(1,g) = xi(0,g)*(1.0 + xi(0,g))*0.50;
447 N(2,g) = (1.0 - xi(0,g))*(1.0 + xi(0,g));
448
449 Nx(0,0,g) = -0.50 + xi(0,g);
450 Nx(0,1,g) = 0.50 + xi(0,g);
451 Nx(0,2,g) = -2.0*xi(0,g);
452 }
453 },
454
455 {ElementType::QUD4, [](const int insd, const int eNoN, const int g, Array<double>& xi, Array<double>& N,
456 Array3<double>& Nx) -> void
457 {
458 double lx = 1.0 - xi(0,g);
459 double ly = 1.0 - xi(1,g);
460 double ux = 1.0 + xi(0,g);
461 double uy = 1.0 + xi(1,g);
462
463 N(0,g) = lx*ly / 4.0;
464 N(1,g) = ux*ly / 4.0;
465 N(2,g) = ux*uy / 4.0;
466 N(3,g) = lx*uy / 4.0;
467
468 Nx(0,0,g) = -ly / 4.0;
469 Nx(1,0,g) = -lx / 4.0;
470 Nx(0,1,g) = ly / 4.0;
471 Nx(1,1,g) = -ux / 4.0;
472 Nx(0,2,g) = uy / 4.0;
473 Nx(1,2,g) = ux / 4.0;
474 Nx(0,3,g) = -uy / 4.0;
475 Nx(1,3,g) = lx / 4.0;
476 }
477 },
478
479 {ElementType::QUD9, [](const int insd, const int eNoN, const int g, Array<double>& xi, Array<double>& N,
480 Array3<double>& Nx) -> void
481 {
482 double lx = 1.0 - xi(0,g);
483 double ly = 1.0 - xi(1,g);
484 double ux = 1.0 + xi(0,g);
485 double uy = 1.0 + xi(1,g);
486 double mx = xi(0,g);
487 double my = xi(1,g);
488
489 N(0,g) = mx*lx*my*ly/4.0;
490 N(1,g) = -mx*ux*my*ly/4.0;
491 N(2,g) = mx*ux*my*uy/4.0;
492 N(3,g) = -mx*lx*my*uy/4.0;
493 N(4,g) = -lx*ux*my*ly*0.50;
494 N(5,g) = mx*ux*ly*uy*0.50;
495 N(6,g) = lx*ux*my*uy*0.50;
496 N(7,g) = -mx*lx*ly*uy*0.50;
497 N(8,g) = lx*ux*ly*uy;
498
499 Nx(0,0,g) = (lx - mx)*my*ly/4.0;
500 Nx(1,0,g) = (ly - my)*mx*lx/4.0;
501 Nx(0,1,g) = -(ux + mx)*my*ly/4.0;
502 Nx(1,1,g) = -(ly - my)*mx*ux/4.0;
503 Nx(0,2,g) = (ux + mx)*my*uy/4.0;
504 Nx(1,2,g) = (uy + my)*mx*ux/4.0;
505 Nx(0,3,g) = -(lx - mx)*my*uy/4.0;
506 Nx(1,3,g) = -(uy + my)*mx*lx/4.0;
507 Nx(0,4,g) = -(lx - ux)*my*ly*0.50;
508 Nx(1,4,g) = -(ly - my)*lx*ux*0.50;
509 Nx(0,5,g) = (ux + mx)*ly*uy*0.50;
510 Nx(1,5,g) = (ly - uy)*mx*ux*0.50;
511 Nx(0,6,g) = (lx - ux)*my*uy*0.50;
512 Nx(1,6,g) = (uy + my)*lx*ux*0.50;
513 Nx(0,7,g) = -(lx - mx)*ly*uy*0.50;
514 Nx(1,7,g) = -(ly - uy)*mx*lx*0.50;
515 Nx(0,8,g) = (lx - ux)*ly*uy;
516 Nx(1,8,g) = (ly - uy)*lx*ux;
517 }
518 },
519
520 {ElementType::TET4, [](const int insd, const int eNoN, const int g, Array<double>& xi, Array<double>& N,
521 Array3<double>& Nx) -> void
522 {
523 //std::cout << "[get_element_shape_data] TET4 " << std::endl;
524
525 N(0,g) = xi(0,g);
526 N(1,g) = xi(1,g);
527 N(2,g) = xi(2,g);
528 N(3,g) = 1.0 - xi(0,g) - xi(1,g) - xi(2,g);
529
530 Nx(0,0,g) = 1.0;
531 Nx(1,0,g) = 0.0;
532 Nx(2,0,g) = 0.0;
533 Nx(0,1,g) = 0.0;
534 Nx(1,1,g) = 1.0;
535 Nx(2,1,g) = 0.0;
536 Nx(0,2,g) = 0.0;
537 Nx(1,2,g) = 0.0;
538 Nx(2,2,g) = 1.0;
539 Nx(0,3,g) = -1.0;
540 Nx(1,3,g) = -1.0;
541 Nx(2,3,g) = -1.0;
542 }
543 },
544
545 {ElementType::TET10, [](const int insd, const int eNoN, const int g, Array<double>& xi, Array<double>& N,
546 Array3<double>& Nx) -> void
547 {
548 double s = 1.0 - xi(0,g) - xi(1,g) - xi(2,g);
549 N(0,g) = xi(0,g)*(2.0*xi(0,g) - 1.0);
550 N(1,g) = xi(1,g)*(2.0*xi(1,g) - 1.0);
551 N(2,g) = xi(2,g)*(2.0*xi(2,g) - 1.0);
552 N(3,g) = s * (2.0*s - 1.0);
553 N(4,g) = 4.0*xi(0,g)*xi(1,g);
554 N(5,g) = 4.0*xi(1,g)*xi(2,g);
555 N(6,g) = 4.0*xi(0,g)*xi(2,g);
556 N(7,g) = 4.0*xi(0,g)*s;
557 N(8,g) = 4.0*xi(1,g)*s;
558 N(9,g) = 4.0*xi(2,g)*s;
559
560 Nx(0,0,g) = 4.0*xi(0,g) - 1.0;
561 Nx(1,0,g) = 0.0;
562 Nx(2,0,g) = 0.0;
563
564 Nx(0,1,g) = 0.0;
565 Nx(1,1,g) = 4.0*xi(1,g) - 1.0;
566 Nx(2,1,g) = 0.0;
567
568 Nx(0,2,g) = 0.0;
569 Nx(1,2,g) = 0.0;
570 Nx(2,2,g) = 4.0*xi(2,g) - 1.0;
571
572 Nx(0,3,g) = 1.0 - 4.0*s;
573 Nx(1,3,g) = 1.0 - 4.0*s;
574 Nx(2,3,g) = 1.0 - 4.0*s;
575
576 Nx(0,4,g) = 4.0*xi(1,g);
577 Nx(1,4,g) = 4.0*xi(0,g);
578 Nx(2,4,g) = 0.0;
579
580 Nx(0,5,g) = 0.0;
581 Nx(1,5,g) = 4.0*xi(2,g);
582 Nx(2,5,g) = 4.0*xi(1,g);
583
584 Nx(0,6,g) = 4.0*xi(2,g);
585 Nx(1,6,g) = 0.0;
586 Nx(2,6,g) = 4.0*xi(0,g);
587
588 Nx(0,7,g) = 4.0*( s - xi(0,g));
589 Nx(1,7,g) = -4.0*xi(0,g);
590 Nx(2,7,g) = -4.0*xi(0,g);
591
592 Nx(0,8,g) = -4.0*xi(1,g);
593 Nx(1,8,g) = 4.0*( s - xi(1,g));
594 Nx(2,8,g) = -4.0*xi(1,g);
595
596 Nx(0,9,g) = -4.0*xi(2,g);
597 Nx(1,9,g) = -4.0*xi(2,g);
598 Nx(2,9,g) = 4.0*( s - xi(2,g));
599 }
600 },
601
602 {ElementType::TRI3, [](const int insd, const int eNoN, const int g, Array<double>& xi, Array<double>& N,
603 Array3<double>& Nx) -> void
604 {
605 //std::cout << "[get_element_shape_data] TRI3 " << std::endl;
606 N(0,g) = xi(0,g);
607 N(1,g) = xi(1,g);
608 N(2,g) = 1.0 - xi(0,g) - xi(1,g);
609
610 Nx(0,0,g) = 1.0;
611 Nx(1,0,g) = 0.0;
612 Nx(0,1,g) = 0.0;
613 Nx(1,1,g) = 1.0;
614 Nx(0,2,g) = -1.0;
615 Nx(1,2,g) = -1.0;
616 }
617 },
618
619 {ElementType::TRI6, [](const int insd, const int eNoN, const int g, Array<double>& xi, Array<double>& N,
620 Array3<double>& Nx) -> void
621 {
622 double s = 1.0 - xi(0,g) - xi(1,g);
623 N(0,g) = xi(0,g) * (2.0*xi(0,g) - 1.0);
624 N(1,g) = xi(1,g) * (2.0*xi(1,g) - 1.0);
625 N(2,g) = s * (2.0*s - 1.0);
626 N(3,g) = 4.0*xi(0,g)*xi(1,g);
627 N(4,g) = 4.0*xi(1,g)*s;
628 N(5,g) = 4.0*xi(0,g)*s;
629
630 Nx(0,0,g) = 4.0*xi(0,g) - 1.0;
631 Nx(1,0,g) = 0.0;
632
633 Nx(0,1,g) = 0.0;
634 Nx(1,1,g) = 4.0*xi(1,g) - 1.0;
635
636 Nx(0,2,g) = 1.0 - 4.0*s;
637 Nx(1,2,g) = 1.0 - 4.0*s;
638
639 Nx(0,3,g) = 4.0*xi(1,g);
640 Nx(1,3,g) = 4.0*xi(0,g);
641
642 Nx(0,4,g) = -4.0*xi(1,g);
643 Nx(1,4,g) = 4.0*( s - xi(1,g) );
644
645 Nx(0,5,g) = 4.0*( s - xi(0,g) );
646 Nx(1,5,g) = -4.0*xi(0,g);
647 }
648 },
649
650 {ElementType::WDG, [](const int insd, const int eNoN, const int g, Array<double>& xi, Array<double>& N,
651 Array3<double>& Nx) -> void
652 {
653 double ux = xi(0,g);
654 double uy = xi(1,g);
655 double uz = 1.0 - ux - uy;
656 double s = (1.0 + xi(2,g))*0.5;
657 double t = (1.0 - xi(2,g))*0.5;
658 N(0,g) = ux*t;
659 N(1,g) = uy*t;
660 N(2,g) = uz*t;
661 N(3,g) = ux*s;
662 N(4,g) = uy*s;
663 N(5,g) = uz*s;
664
665 Nx(0,0,g) = t;
666 Nx(1,0,g) = 0.0;
667 Nx(2,0,g) = -ux*0.50;
668
669 Nx(0,1,g) = 0.0;
670 Nx(1,1,g) = t;
671 Nx(2,1,g) = -uy*0.50;
672
673 Nx(0,2,g) = -t;
674 Nx(1,2,g) = -t;
675 Nx(2,2,g) = -uz*0.50;
676
677 Nx(0,3,g) = s;
678 Nx(1,3,g) = 0.0;
679 Nx(2,3,g) = ux*0.50;
680
681 Nx(0,4,g) = 0.0;
682 Nx(1,4,g) = s;
683 Nx(2,4,g) = uy*0.50;
684
685 Nx(0,5,g) = -s;
686 Nx(1,5,g) = -s;
687 Nx(2,5,g) = uz*0.50;
688 }
689 },
690
691
692
693};
694
695
696//------------------------
697// set_element_shape_data
698//------------------------
699// Replicates 'SUBROUTINE GETGNN(insd, eType, eNoN, xi, N, Nxi)' defined in NN.f.
700//
701using SetElementShapeMapType = std::map<ElementType, std::function<void(int, mshType&)>>;
702
703SetElementShapeMapType set_element_shape_data = {
704
705 {ElementType::HEX8, [](int g, mshType& mesh) -> void {
706 auto& xi = mesh.xi;
707 double lx = 1.0 - xi(0,g);
708 double ly = 1.0 - xi(1,g);
709 double lz = 1.0 - xi(2,g);
710 double ux = 1.0 + xi(0,g);
711 double uy = 1.0 + xi(1,g);
712 double uz = 1.0 + xi(2,g);
713
714 auto& N = mesh.N;
715 N(0,g) = lx*ly*lz/8.0;
716 N(1,g) = ux*ly*lz/8.0;
717 N(2,g) = ux*uy*lz/8.0;
718 N(3,g) = lx*uy*lz/8.0;
719 N(4,g) = lx*ly*uz/8.0;
720 N(5,g) = ux*ly*uz/8.0;
721 N(6,g) = ux*uy*uz/8.0;
722 N(7,g) = lx*uy*uz/8.0;
723
724 auto& Nx = mesh.Nx;
725 Nx(0,0,g) = -ly*lz/8.0;
726 Nx(1,0,g) = -lx*lz/8.0;
727 Nx(2,0,g) = -lx*ly/8.0;
728
729 Nx(0,1,g) = ly*lz/8.0;
730 Nx(1,1,g) = -ux*lz/8.0;
731 Nx(2,1,g) = -ux*ly/8.0;
732
733 Nx(0,2,g) = uy*lz/8.0;
734 Nx(1,2,g) = ux*lz/8.0;
735 Nx(2,2,g) = -ux*uy/8.0;
736
737 Nx(0,3,g) = -uy*lz/8.0;
738 Nx(1,3,g) = lx*lz/8.0;
739 Nx(2,3,g) = -lx*uy/8.0;
740
741 Nx(0,4,g) = -ly*uz/8.0;
742 Nx(1,4,g) = -lx*uz/8.0;
743 Nx(2,4,g) = lx*ly/8.0;
744
745 Nx(0,5,g) = ly*uz/8.0;
746 Nx(1,5,g) = -ux*uz/8.0;
747 Nx(2,5,g) = ux*ly/8.0;
748
749 Nx(0,6,g) = uy*uz/8.0;
750 Nx(1,6,g) = ux*uz/8.0;
751 Nx(2,6,g) = ux*uy/8.0;
752
753 Nx(0,7,g) = -uy*uz/8.0;
754 Nx(1,7,g) = lx*uz/8.0;
755 Nx(2,7,g) = lx*uy/8.0;
756 }
757 },
758
759 {ElementType::HEX20, [](int g, mshType& mesh) -> void {
760
761 auto& xi = mesh.xi;
762 double lx = 1.0 - xi(0,g);
763 double ly = 1.0 - xi(1,g);
764 double lz = 1.0 - xi(2,g);
765 double ux = 1.0 + xi(0,g);
766 double uy = 1.0 + xi(1,g);
767 double uz = 1.0 + xi(2,g);
768
769 double mx = lx*ux;
770 double my = ly*uy;
771 double mz = lz*uz;
772
773 auto& N = mesh.N;
774 N(0, g) = lx*ly*lz*(lx+ly+lz-5.0)/8.0;
775 N(1, g) = ux*ly*lz*(ux+ly+lz-5.0)/8.0;
776 N(2, g) = ux*uy*lz*(ux+uy+lz-5.0)/8.0;
777 N(3, g) = lx*uy*lz*(lx+uy+lz-5.0)/8.0;
778 N(4, g) = lx*ly*uz*(lx+ly+uz-5.0)/8.0;
779 N(5, g) = ux*ly*uz*(ux+ly+uz-5.0)/8.0;
780 N(6, g) = ux*uy*uz*(ux+uy+uz-5.0)/8.0;
781 N(7, g) = lx*uy*uz*(lx+uy+uz-5.0)/8.0;
782 N(8, g) = mx*ly*lz/4.0;
783 N(9, g) = ux*my*lz/4.0;
784 N(10, g) = mx*uy*lz/4.0;
785 N(11, g) = lx*my*lz/4.0;
786 N(12, g) = mx*ly*uz/4.0;
787 N(13, g) = ux*my*uz/4.0;
788 N(14, g) = mx*uy*uz/4.0;
789 N(15, g) = lx*my*uz/4.0;
790 N(16, g) = lx*ly*mz/4.0;
791 N(17, g) = ux*ly*mz/4.0;
792 N(18, g) = ux*uy*mz/4.0;
793 N(19, g) = lx*uy*mz/4.0;
794
795 // N(1) = lx*ly*lz*(lx+ly+lz-5.0)/8.0;
796 auto& Nx = mesh.Nx;
797 int n = 0;
798 Nx(0,n,g) = -ly*lz*(lx+ly+lz-5.0+lx)/8.0;
799 Nx(1,n,g) = -lx*lz*(lx+ly+lz-5.0+ly)/8.0;
800 Nx(2,n,g) = -lx*ly*(lx+ly+lz-5.0+lz)/8.0;
801
802//c N(n,g) = ux*ly*lz*(ux+ly+lz-5.0)/8.0;
803 n += 1;
804 Nx(0,n,g) = ly*lz*(ux+ly+lz-5.0+ux)/8.0;
805 Nx(1,n,g) = -ux*lz*(ux+ly+lz-5.0+ly)/8.0;
806 Nx(2,n,g) = -ux*ly*(ux+ly+lz-5.0+lz)/8.0;
807
808//c N(n,g) = ux*uy*lz*(ux+uy+lz-5.0)/8.0
809 n += 1;
810 Nx(0,n,g) = uy*lz*(ux+uy+lz-5.0+ux)/8.0;
811 Nx(1,n,g) = ux*lz*(ux+uy+lz-5.0+uy)/8.0;
812 Nx(2,n,g) = -ux*uy*(ux+uy+lz-5.0+lz)/8.0;
813
814//c N(n,g) = lx*uy*lz*(lx+uy+lz-5.0)/8.0
815 n += 1;
816 Nx(0,n,g) = -uy*lz*(lx+uy+lz-5.0+lx)/8.0;
817 Nx(1,n,g) = lx*lz*(lx+uy+lz-5.0+uy)/8.0;
818 Nx(2,n,g) = -lx*uy*(lx+uy+lz-5.0+lz)/8.0;
819
820//c N(n,g) = lx*ly*uz*(lx+ly+uz-5.0)/8.0
821 n += 1;
822 Nx(0,n,g) = -ly*uz*(lx+ly+uz-5.0+lx)/8.0;
823 Nx(1,n,g) = -lx*uz*(lx+ly+uz-5.0+ly)/8.0;
824 Nx(2,n,g) = lx*ly*(lx+ly+uz-5.0+uz)/8.0;
825
826//c N(n,g) = ux*ly*uz*(ux+ly+uz-5.0)/8.0
827 n += 1;
828 Nx(0,n,g) = ly*uz*(ux+ly+uz-5.0+ux)/8.0;
829 Nx(1,n,g) = -ux*uz*(ux+ly+uz-5.0+ly)/8.0;
830 Nx(2,n,g) = ux*ly*(ux+ly+uz-5.0+uz)/8.0;
831
832//c N(n,g) = ux*uy*uz*(ux+uy+uz-5.0)/8.0
833 n += 1;
834 Nx(0,n,g) = uy*uz*(ux+uy+uz-5.0+ux)/8.0;
835 Nx(1,n,g) = ux*uz*(ux+uy+uz-5.0+uy)/8.0;
836 Nx(2,n,g) = ux*uy*(ux+uy+uz-5.0+uz)/8.0;
837
838//c N(n,g) = lx*uy*uz*(lx+uy+uz-5.0)/8.0
839 n += 1;
840 Nx(0,n,g) = -uy*uz*(lx+uy+uz-5.0+lx)/8.0;
841 Nx(1,n,g) = lx*uz*(lx+uy+uz-5.0+uy)/8.0;
842 Nx(2,n,g) = lx*uy*(lx+uy+uz-5.0+uz)/8.0;
843
844//c N(n,g) = mx*ly*lz/4.0
845 n += 1;
846 Nx(0,n,g) = (lx - ux)*ly*lz/4.0;
847 Nx(1,n,g) = -mx*lz/4.0;
848 Nx(2,n,g) = -mx*ly/4.0;
849
850//c N(0n,g) = ux*my*lz/4.0
851 n += 1;
852 Nx(0,n,g) = my*lz/4.0;
853 Nx(1,n,g) = (ly - uy)*ux*lz/4.0;
854 Nx(2,n,g) = -ux*my/4.0;
855
856//c N(0n,g) = mx*uy*lz/4.0
857 n += 1;
858 Nx(0,n,g) = (lx - ux)*uy*lz/4.0;
859 Nx(1,n,g) = mx*lz/4.0;
860 Nx(2,n,g) = -mx*uy/4.0;
861
862//c N(0n,g) = lx*my*lz/4.0
863 n += 1;
864 Nx(0,n,g) = -my*lz/4.0;
865 Nx(1,n,g) = (ly - uy)*lx*lz/4.0;
866 Nx(2,n,g) = -lx*my/4.0;
867
868//c N(0n,g) = mx*ly*uz/4.0
869 n += 1;
870 Nx(0,n,g) = (lx - ux)*ly*uz/4.0;
871 Nx(1,n,g) = -mx*uz/4.0;
872 Nx(2,n,g) = mx*ly/4.0;
873
874//c N(0n,g) = ux*my*uz/4.0
875 n += 1;
876 Nx(0,n,g) = my*uz/4.0;
877 Nx(1,n,g) = (ly - uy)*ux*uz/4.0;
878 Nx(2,n,g) = ux*my/4.0;
879
880//c N(0n,g) = mx*uy*uz/4.0
881 n += 1;
882 Nx(0,n,g) = (lx - ux)*uy*uz/4.0;
883 Nx(1,n,g) = mx*uz/4.0;
884 Nx(2,n,g) = mx*uy/4.0;
885
886//c N(0n,g) = lx*my*uz/4.0
887 n += 1;
888 Nx(0,n,g) = -my*uz/4.0;
889 Nx(1,n,g) = (ly - uy)*lx*uz/4.0;
890 Nx(2,n,g) = lx*my/4.0;
891
892//c N(0n,g) = lx*ly*mz/4.0
893 n += 1;
894 Nx(0,n,g) = -ly*mz/4.0;
895 Nx(1,n,g) = -lx*mz/4.0;
896 Nx(2,n,g) = (lz - uz)*lx*ly/4.0;
897
898//c N(0n,g) = ux*ly*mz/4.0
899 n += 1;
900 Nx(0,n,g) = ly*mz/4.0;
901 Nx(1,n,g) = -ux*mz/4.0;
902 Nx(2,n,g) = (lz - uz)*ux*ly/4.0;
903
904//c N(0n,g) = ux*uy*mz/4.0
905 n += 1;
906 Nx(0,n,g) = uy*mz/4.0;
907 Nx(1,n,g) = ux*mz/4.0;
908 Nx(2,n,g) = (lz - uz)*ux*uy/4.0;
909
910//c N(n,g) = lx*uy*mz/4.0
911 n += 1;
912 Nx(0,n,g) = -uy*mz/4.0;
913 Nx(1,n,g) = lx*mz/4.0;
914 Nx(2,n,g) = (lz - uz)*lx*uy/4.0;
915 }
916 },
917
918 {ElementType::HEX27, [](int g, mshType& mesh) -> void {
919
920 auto& xi = mesh.xi;
921 double lx = 1.0 - xi(0,g);
922 double ly = 1.0 - xi(1,g);
923 double lz = 1.0 - xi(2,g);
924 double ux = 1.0 + xi(0,g);
925 double uy = 1.0 + xi(1,g);
926 double uz = 1.0 + xi(2,g);
927
928 double mx = xi(0,g);
929 double my = xi(1,g);
930 double mz = xi(2,g);
931
932 auto& N = mesh.N;
933 N(0,g) = -mx*lx*my*ly*mz*lz/8.0;
934 N(1,g) = mx*ux*my*ly*mz*lz/8.0;
935 N(2,g) = -mx*ux*my*uy*mz*lz/8.0;
936 N(3,g) = mx*lx*my*uy*mz*lz/8.0;
937 N(4,g) = mx*lx*my*ly*mz*uz/8.0;
938 N(5,g) = -mx*ux*my*ly*mz*uz/8.0;
939 N(6,g) = mx*ux*my*uy*mz*uz/8.0;
940 N(7,g) = -mx*lx*my*uy*mz*uz/8.0;
941 N(8,g) = lx*ux*my*ly*mz*lz/4.0;
942 N(9,g) = -mx*ux*ly*uy*mz*lz/4.0;
943 N(10,g) = -lx*ux*my*uy*mz*lz/4.0;
944 N(11,g) = mx*lx*ly*uy*mz*lz/4.0;
945 N(12,g) = -lx*ux*my*ly*mz*uz/4.0;
946 N(13,g) = mx*ux*ly*uy*mz*uz/4.0;
947 N(14,g) = lx*ux*my*uy*mz*uz/4.0;
948 N(15,g) = -mx*lx*ly*uy*mz*uz/4.0;
949 N(16,g) = mx*lx*my*ly*lz*uz/4.0;
950 N(17,g) = -mx*ux*my*ly*lz*uz/4.0;
951 N(18,g) = mx*ux*my*uy*lz*uz/4.0;
952 N(19,g) = -mx*lx*my*uy*lz*uz/4.0;
953
954 N(20,g) = -mx*lx*ly*uy*lz*uz/2.0;
955 N(21,g) = mx*ux*ly*uy*lz*uz/2.0;
956 N(22,g) = -lx*ux*my*ly*lz*uz/2.0;
957 N(23,g) = lx*ux*my*uy*lz*uz/2.0;
958 N(24,g) = -lx*ux*ly*uy*mz*lz/2.0;
959 N(25,g) = lx*ux*ly*uy*mz*uz/2.0;
960
961 N(26,g) = lx*ux*ly*uy*lz*uz;
962
963 auto& Nxi = mesh.Nx;
964 int n = 0;
965 Nxi(0,n,g) = -(lx - mx)*my*ly*mz*lz/8.0;
966 Nxi(1,n,g) = -(ly - my)*mx*lx*mz*lz/8.0;
967 Nxi(2,n,g) = -(lz - mz)*mx*lx*my*ly/8.0;
968
969 n += 1;
970 Nxi(0,n,g) = (mx + ux)*my*ly*mz*lz/8.0;
971 Nxi(1,n,g) = (ly - my)*mx*ux*mz*lz/8.0;
972 Nxi(2,n,g) = (lz - mz)*mx*ux*my*ly/8.0;
973
974 n += 1;
975 Nxi(0,n,g) = -(mx + ux)*my*uy*mz*lz/8.0;
976 Nxi(1,n,g) = -(my + uy)*mx*ux*mz*lz/8.0;
977 Nxi(2,n,g) = -(lz - mz)*mx*ux*my*uy/8.0;
978
979 n += 1;
980 Nxi(0,n,g) = (lx - mx)*my*uy*mz*lz/8.0;
981 Nxi(1,n,g) = (my + uy)*mx*lx*mz*lz/8.0;
982 Nxi(2,n,g) = (lz - mz)*mx*lx*my*uy/8.0;
983
984 n += 1;
985 Nxi(0,n,g) = (lx - mx)*my*ly*mz*uz/8.0;
986 Nxi(1,n,g) = (ly - my)*mx*lx*mz*uz/8.0;
987 Nxi(2,n,g) = (mz + uz)*mx*lx*my*ly/8.0;
988
989 n += 1;
990 Nxi(0,n,g) = -(mx + ux)*my*ly*mz*uz/8.0;
991 Nxi(1,n,g) = -(ly - my)*mx*ux*mz*uz/8.0;
992 Nxi(2,n,g) = -(mz + uz)*mx*ux*my*ly/8.0;
993
994 n += 1;
995 Nxi(0,n,g) = (mx + ux)*my*uy*mz*uz/8.0;
996 Nxi(1,n,g) = (my + uy)*mx*ux*mz*uz/8.0;
997 Nxi(2,n,g) = (mz + uz)*mx*ux*my*uy/8.0;
998
999 n += 1;
1000 Nxi(0,n,g) = -(lx - mx)*my*uy*mz*uz/8.0;
1001 Nxi(1,n,g) = -(my + uy)*mx*lx*mz*uz/8.0;
1002 Nxi(2,n,g) = -(mz + uz)*mx*lx*my*uy/8.0;
1003
1004 n += 1;
1005 Nxi(0,n,g) = (lx - ux)*my*ly*mz*lz/4.0;
1006 Nxi(1,n,g) = (ly - my)*lx*ux*mz*lz/4.0;
1007 Nxi(2,n,g) = (lz - mz)*lx*ux*my*ly/4.0;
1008
1009 n += 1;
1010 Nxi(0,n,g) = -(mx + ux)*ly*uy*mz*lz/4.0;
1011 Nxi(1,n,g) = -(ly - uy)*mx*ux*mz*lz/4.0;
1012 Nxi(2,n,g) = -(lz - mz)*mx*ux*ly*uy/4.0;
1013
1014 n += 1;
1015 Nxi(0,n,g) = -(lx - ux)*my*uy*mz*lz/4.0;
1016 Nxi(1,n,g) = -(my + uy)*lx*ux*mz*lz/4.0;
1017 Nxi(2,n,g) = -(lz - mz)*lx*ux*my*uy/4.0;
1018
1019 n += 1;
1020 Nxi(0,n,g) = (lx - mx)*ly*uy*mz*lz/4.0;
1021 Nxi(1,n,g) = (ly - uy)*mx*lx*mz*lz/4.0;
1022 Nxi(2,n,g) = (lz - mz)*mx*lx*ly*uy/4.0;
1023
1024 n += 1;
1025 Nxi(0,n,g) = -(lx - ux)*my*ly*mz*uz/4.0;
1026 Nxi(1,n,g) = -(ly - my)*lx*ux*mz*uz/4.0;
1027 Nxi(2,n,g) = -(mz + uz)*lx*ux*my*ly/4.0;
1028
1029 n += 1;
1030 Nxi(0,n,g) = (mx + ux)*ly*uy*mz*uz/4.0;
1031 Nxi(1,n,g) = (ly - uy)*mx*ux*mz*uz/4.0;
1032 Nxi(2,n,g) = (mz + uz)*mx*ux*ly*uy/4.0;
1033
1034 n += 1;
1035 Nxi(0,n,g) = (lx - ux)*my*uy*mz*uz/4.0;
1036 Nxi(1,n,g) = (my + uy)*lx*ux*mz*uz/4.0;
1037 Nxi(2,n,g) = (mz + uz)*lx*ux*my*uy/4.0;
1038
1039 n += 1;
1040 Nxi(0,n,g) = -(lx - mx)*ly*uy*mz*uz/4.0;
1041 Nxi(1,n,g) = -(ly - uy)*mx*lx*mz*uz/4.0;
1042 Nxi(2,n,g) = -(mz + uz)*mx*lx*ly*uy/4.0;
1043
1044 n += 1;
1045 Nxi(0,n,g) = (lx - mx)*my*ly*lz*uz/4.0;
1046 Nxi(1,n,g) = (ly - my)*mx*lx*lz*uz/4.0;
1047 Nxi(2,n,g) = (lz - uz)*mx*lx*my*ly/4.0;
1048
1049 n += 1;
1050 Nxi(0,n,g) = -(mx + ux)*my*ly*lz*uz/4.0;
1051 Nxi(1,n,g) = -(ly - my)*mx*ux*lz*uz/4.0;
1052 Nxi(2,n,g) = -(lz - uz)*mx*ux*my*ly/4.0;
1053
1054 n += 1;
1055 Nxi(0,n,g) = (mx + ux)*my*uy*lz*uz/4.0;
1056 Nxi(1,n,g) = (my + uy)*mx*ux*lz*uz/4.0;
1057 Nxi(2,n,g) = (lz - uz)*mx*ux*my*uy/4.0;
1058
1059 n += 1;
1060 Nxi(0,n,g) = -(lx - mx)*my*uy*lz*uz/4.0;
1061 Nxi(1,n,g) = -(my + uy)*mx*lx*lz*uz/4.0;
1062 Nxi(2,n,g) = -(lz - uz)*mx*lx*my*uy/4.0;
1063
1064 n += 1;
1065 Nxi(0,n,g) = -(lx - mx)*ly*uy*lz*uz/2.0;
1066 Nxi(1,n,g) = -(ly - uy)*mx*lx*lz*uz/2.0;
1067 Nxi(2,n,g) = -(lz - uz)*mx*lx*ly*uy/2.0;
1068
1069 n += 1;
1070 Nxi(0,n,g) = (mx + ux)*ly*uy*lz*uz/2.0;
1071 Nxi(1,n,g) = (ly - uy)*mx*ux*lz*uz/2.0;
1072 Nxi(2,n,g) = (lz - uz)*mx*ux*ly*uy/2.0;
1073
1074 n += 1;
1075 Nxi(0,n,g) = -(lx - ux)*my*ly*lz*uz/2.0;
1076 Nxi(1,n,g) = -(ly - my)*lx*ux*lz*uz/2.0;
1077 Nxi(2,n,g) = -(lz - uz)*lx*ux*my*ly/2.0;
1078
1079 n += 1;
1080 Nxi(0,n,g) = (lx - ux)*my*uy*lz*uz/2.0;
1081 Nxi(1,n,g) = (my + uy)*lx*ux*lz*uz/2.0;
1082 Nxi(2,n,g) = (lz - uz)*lx*ux*my*uy/2.0;
1083
1084 n += 1;
1085 Nxi(0,n,g) = -(lx - ux)*ly*uy*mz*lz/2.0;
1086 Nxi(1,n,g) = -(ly - uy)*lx*ux*mz*lz/2.0;
1087 Nxi(2,n,g) = -(lz - mz)*lx*ux*ly*uy/2.0;
1088
1089 n += 1;
1090 Nxi(0,n,g) = (lx - ux)*ly*uy*mz*uz/2.0;
1091 Nxi(1,n,g) = (ly - uy)*lx*ux*mz*uz/2.0;
1092 Nxi(2,n,g) = (mz + uz)*lx*ux*ly*uy/2.0;
1093
1094 n += 1;
1095 Nxi(0,n,g) = (lx - ux)*ly*uy*lz*uz;
1096 Nxi(1,n,g) = (ly - uy)*lx*ux*lz*uz;
1097 Nxi(2,n,g) = (lz - uz)*lx*ux*ly*uy;
1098 }
1099 },
1100
1101 {ElementType::LIN1, [](int g, mshType& mesh) -> void {
1102 //std::cout << "[set_element_shape_data] **************************" << std::endl;
1103 //std::cout << "[set_element_shape_data] ERROR: LIN1 not supported." << std::endl;
1104 //std::cout << "[set_element_shape_data] **************************" << std::endl;
1105 auto& xi = mesh.xi;
1106 auto& N = mesh.N;
1107 N(0,g) = (1.0 - xi(0,g))*0.5;
1108 N(1,g) = (1.0 + xi(0,g))*0.5;
1109
1110 auto& Nx = mesh.Nx;
1111 Nx(0,0,g) = -0.5;
1112 Nx(0,1,g) = 0.5;
1113 }
1114 },
1115
1116 {ElementType::LIN2, [](int g, mshType& mesh) -> void {
1117 auto& xi = mesh.xi;
1118 auto& N = mesh.N;
1119 N(0,g) = -xi(0,g)*(1.0 - xi(0,g))*0.50;
1120 N(1,g) = xi(0,g)*(1.0 + xi(0,g))*0.50;
1121 N(2,g) = (1.0 - xi(0,g))*(1.0 + xi(0,g));
1122
1123 auto& Nx = mesh.Nx;
1124 Nx(0,0,g) = -0.50 + xi(0,g);
1125 Nx(0,1,g) = 0.50 + xi(0,g);
1126 Nx(0,2,g) = -2.0*xi(0,g);
1127 }
1128 },
1129
1130 {ElementType::QUD4, [](int g, mshType& mesh) -> void {
1131 auto& xi = mesh.xi;
1132 double lx = 1.0 - xi(0,g);
1133 double ly = 1.0 - xi(1,g);
1134 double ux = 1.0 + xi(0,g);
1135 double uy = 1.0 + xi(1,g);
1136
1137 auto& N = mesh.N;
1138 N(0,g) = lx*ly / 4.0;
1139 N(1,g) = ux*ly / 4.0;
1140 N(2,g) = ux*uy / 4.0;
1141 N(3,g) = lx*uy / 4.0;
1142
1143 auto& Nx = mesh.Nx;
1144 Nx(0,0,g) = -ly / 4.0;
1145 Nx(1,0,g) = -lx / 4.0;
1146 Nx(0,1,g) = ly / 4.0;
1147 Nx(1,1,g) = -ux / 4.0;
1148 Nx(0,2,g) = uy / 4.0;
1149 Nx(1,2,g) = ux / 4.0;
1150 Nx(0,3,g) = -uy / 4.0;
1151 Nx(1,3,g) = lx / 4.0;
1152 }
1153 },
1154
1155 {ElementType::QUD9, [](int g, mshType& mesh) -> void {
1156 auto& xi = mesh.xi;
1157 double lx = 1.0 - xi(0,g);
1158 double ly = 1.0 - xi(1,g);
1159 double ux = 1.0 + xi(0,g);
1160 double uy = 1.0 + xi(1,g);
1161 double mx = xi(0,g);
1162 double my = xi(1,g);
1163
1164 auto& N = mesh.N;
1165 N(0,g) = mx*lx*my*ly/4.0;
1166 N(1,g) = -mx*ux*my*ly/4.0;
1167 N(2,g) = mx*ux*my*uy/4.0;
1168 N(3,g) = -mx*lx*my*uy/4.0;
1169 N(4,g) = -lx*ux*my*ly*0.50;
1170 N(5,g) = mx*ux*ly*uy*0.50;
1171 N(6,g) = lx*ux*my*uy*0.50;
1172 N(7,g) = -mx*lx*ly*uy*0.50;
1173 N(8,g) = lx*ux*ly*uy;
1174
1175 auto& Nx = mesh.Nx;
1176 Nx(0,0,g) = (lx - mx)*my*ly/4.0;
1177 Nx(1,0,g) = (ly - my)*mx*lx/4.0;
1178
1179 Nx(0,1,g) = -(ux + mx)*my*ly/4.0;
1180 Nx(1,1,g) = -(ly - my)*mx*ux/4.0;
1181
1182 Nx(0,2,g) = (ux + mx)*my*uy/4.0;
1183 Nx(1,2,g) = (uy + my)*mx*ux/4.0;
1184
1185 Nx(0,3,g) = -(lx - mx)*my*uy/4.0;
1186 Nx(1,3,g) = -(uy + my)*mx*lx/4.0;
1187
1188 Nx(0,4,g) = -(lx - ux)*my*ly*0.50;
1189 Nx(1,4,g) = -(ly - my)*lx*ux*0.50;
1190
1191 Nx(0,5,g) = (ux + mx)*ly*uy*0.50;
1192 Nx(1,5,g) = (ly - uy)*mx*ux*0.50;
1193
1194 Nx(0,6,g) = (lx - ux)*my*uy*0.50;
1195 Nx(1,6,g) = (uy + my)*lx*ux*0.50;
1196
1197 Nx(0,7,g) = -(lx - mx)*ly*uy*0.50;
1198 Nx(1,7,g) = -(ly - uy)*mx*lx*0.50;
1199
1200 Nx(0,8,g) = (lx - ux)*ly*uy;
1201 Nx(1,8,g) = (ly - uy)*lx*ux;
1202 }
1203 },
1204
1205 {ElementType::TET4, [](int g, mshType& mesh) -> void {
1206 auto& xi = mesh.xi;
1207 auto& N = mesh.N;
1208 N(0,g) = xi(0,g);
1209 N(1,g) = xi(1,g);
1210 N(2,g) = xi(2,g);
1211 N(3,g) = 1.0 - xi(0,g) - xi(1,g) - xi(2,g);
1212
1213 auto& Nx = mesh.Nx;
1214 Nx(0,0,g) = 1.0;
1215 Nx(1,0,g) = 0.0;
1216 Nx(2,0,g) = 0.0;
1217 Nx(0,1,g) = 0.0;
1218 Nx(1,1,g) = 1.0;
1219 Nx(2,1,g) = 0.0;
1220 Nx(0,2,g) = 0.0;
1221 Nx(1,2,g) = 0.0;
1222 Nx(2,2,g) = 1.0;
1223 Nx(0,3,g) = -1.0;
1224 Nx(1,3,g) = -1.0;
1225 Nx(2,3,g) = -1.0;
1226 }
1227 },
1228
1229 {ElementType::TET10, [](int g, mshType& mesh) -> void {
1230 auto& xi = mesh.xi;
1231 auto& N = mesh.N;
1232 double s = 1.0 - xi(0,g) - xi(1,g) - xi(2,g);
1233 N(0,g) = xi(0,g)*(2.0*xi(0,g) - 1.0);
1234 N(1,g) = xi(1,g)*(2.0*xi(1,g) - 1.0);
1235 N(2,g) = xi(2,g)*(2.0*xi(2,g) - 1.0);
1236 N(3,g) = s *(2.0*s - 1.0);
1237 N(4,g) = 4.0*xi(0,g)*xi(1,g);
1238 N(5,g) = 4.0*xi(1,g)*xi(2,g);
1239 N(6,g) = 4.0*xi(0,g)*xi(2,g);
1240 N(7,g) = 4.0*xi(0,g)*s;
1241 N(8,g) = 4.0*xi(1,g)*s;
1242 N(9,g) = 4.0*xi(2,g)*s;
1243
1244 auto& Nx = mesh.Nx;
1245 Nx(0,0,g) = 4.0*xi(0,g) - 1.0;
1246 Nx(1,0,g) = 0.0;
1247 Nx(2,0,g) = 0.0;
1248
1249 Nx(0,1,g) = 0.0;
1250 Nx(1,1,g) = 4.0*xi(1,g) - 1.0;
1251 Nx(2,1,g) = 0.0;
1252
1253 Nx(0,2,g) = 0.0;
1254 Nx(1,2,g) = 0.0;
1255 Nx(2,2,g) = 4.0*xi(2,g) - 1.0;
1256
1257 Nx(0,3,g) = 1.0 - 4.0*s;
1258 Nx(1,3,g) = 1.0 - 4.0*s;
1259 Nx(2,3,g) = 1.0 - 4.0*s;
1260
1261 Nx(0,4,g) = 4.0*xi(1,g);
1262 Nx(1,4,g) = 4.0*xi(0,g);
1263 Nx(2,4,g) = 0.0;
1264
1265 Nx(0,5,g) = 0.0;
1266 Nx(1,5,g) = 4.0*xi(2,g);
1267 Nx(2,5,g) = 4.0*xi(1,g);
1268
1269 Nx(0,6,g) = 4.0*xi(2,g);
1270 Nx(1,6,g) = 0.0;
1271 Nx(2,6,g) = 4.0*xi(0,g);
1272
1273 Nx(0,7,g) = 4.0*( s - xi(0,g));
1274 Nx(1,7,g) = -4.0*xi(0,g);
1275 Nx(2,7,g) = -4.0*xi(0,g);
1276
1277 Nx(0,8,g) = -4.0*xi(1,g);
1278 Nx(1,8,g) = 4.0*( s - xi(1,g));
1279 Nx(2,8,g) = -4.0*xi(1,g);
1280
1281 Nx(0,9,g) = -4.0*xi(2,g);
1282 Nx(1,9,g) = -4.0*xi(2,g);
1283 Nx(2,9,g) = 4.0*( s - xi(2,g));
1284 }
1285 },
1286
1287 {ElementType::TRI3, [](int g, mshType& mesh) -> void {
1288 auto& xi = mesh.xi;
1289 auto& N = mesh.N;
1290 N(0,g) = xi(0,g);
1291 N(1,g) = xi(1,g);
1292 N(2,g) = 1.0 - xi(0,g) - xi(1,g);
1293
1294 auto& Nxi = mesh.Nx;
1295 Nxi(0,0,g) = 1.0;
1296 Nxi(1,0,g) = 0.0;
1297 Nxi(0,1,g) = 0.0;
1298 Nxi(1,1,g) = 1.0;
1299 Nxi(0,2,g) = -1.0;
1300 Nxi(1,2,g) = -1.0;
1301 }
1302 },
1303
1304 {ElementType::TRI6, [](int g, mshType& mesh) -> void {
1305 auto& xi = mesh.xi;
1306 auto& N = mesh.N;
1307
1308 double s = 1.0 - xi(0,g) - xi(1,g);
1309 N(0,g) = xi(0,g)*( 2.0*xi(0,g) - 1.0 );
1310 N(1,g) = xi(1,g)*( 2.0*xi(1,g) - 1.0 );
1311 N(2,g) = s *( 2.0*s - 1.0 );
1312 N(3,g) = 4.0*xi(0,g)*xi(1,g);
1313 N(4,g) = 4.0*xi(1,g)*s;
1314 N(5,g) = 4.0*xi(0,g)*s;
1315
1316 auto& Nxi = mesh.Nx;
1317 Nxi(0,0,g) = 4.0*xi(0,g) - 1.0;
1318 Nxi(1,0,g) = 0.0;
1319 Nxi(0,1,g) = 0.0;
1320 Nxi(1,1,g) = 4.0*xi(1,g) - 1.0;
1321 Nxi(0,2,g) = 1.0 - 4.0*s;
1322 Nxi(1,2,g) = 1.0 - 4.0*s;
1323 Nxi(0,3,g) = 4.0*xi(1,g);
1324 Nxi(1,3,g) = 4.0*xi(0,g);
1325 Nxi(0,4,g) = -4.0*xi(1,g);
1326 Nxi(1,4,g) = 4.0*( s - xi(1,g) );
1327 Nxi(0,5,g) = 4.0*( s - xi(0,g) );
1328 Nxi(1,5,g) = -4.0*xi(0,g);
1329 }
1330 },
1331
1332 {ElementType::WDG, [](int g, mshType& mesh) -> void
1333 {
1334 auto& xi = mesh.xi;
1335 auto& N = mesh.N;
1336 double ux = xi(0,g);
1337 double uy = xi(1,g);
1338 double uz = 1.0 - ux - uy;
1339 double s = (1.0 + xi(2,g))*0.5;
1340 double t = (1.0 - xi(2,g))*0.5;
1341 N(0,g) = ux*t;
1342 N(1,g) = uy*t;
1343 N(2,g) = uz*t;
1344 N(3,g) = ux*s;
1345 N(4,g) = uy*s;
1346 N(5,g) = uz*s;
1347
1348 auto& Nxi = mesh.Nx;
1349 Nxi(0,0,g) = t;
1350 Nxi(1,0,g) = 0.0;
1351 Nxi(2,0,g) = -ux*0.50;
1352
1353 Nxi(0,1,g) = 0.0;
1354 Nxi(1,1,g) = t;
1355 Nxi(2,1,g) = -uy*0.50;
1356
1357 Nxi(0,2,g) = -t;
1358 Nxi(1,2,g) = -t;
1359 Nxi(2,2,g) = -uz*0.50;
1360
1361 Nxi(0,3,g) = s;
1362 Nxi(1,3,g) = 0.0;
1363 Nxi(2,3,g) = ux*0.50;
1364
1365 Nxi(0,4,g) = 0.0;
1366 Nxi(1,4,g) = s;
1367 Nxi(2,4,g) = uy*0.50;
1368
1369 Nxi(0,5,g) = -s;
1370 Nxi(1,5,g) = -s;
1371 Nxi(2,5,g) = uz*0.50;
1372 }
1373 },
1374
1375};
1376
1377//---------------------
1378// set_face_shape_data
1379//---------------------
1380// Define a map type used to face element shape function data.
1381//
1382// This reproduces 'SUBROUTINE GETGNN(insd, eType, eNoN, xi, N, Nxi)' in NN.f.
1383//
1384using SetFaceShapeMapType = std::map<ElementType, std::function<void(int, faceType&)>>;
1385
1386SetFaceShapeMapType set_face_shape_data = {
1387
1388 {ElementType::PNT, [](int g, faceType& face) -> void
1389 {
1390 face.N(0,g) = 1.0;
1391 }
1392 },
1393
1394 {ElementType::QUD8, [](int g, faceType& face) -> void
1395 {
1396 auto& xi = face.xi;
1397 double lx = 1.0 - xi(0,g);
1398 double ly = 1.0 - xi(1,g);
1399 double ux = 1.0 + xi(0,g);
1400 double uy = 1.0 + xi(1,g);
1401 double mx = lx*ux;
1402 double my = ly*uy;
1403
1404 auto& N = face.N;
1405 N(0,g) = lx*ly*(lx+ly-3.0)/4.0;
1406 N(1,g) = ux*ly*(ux+ly-3.0)/4.0;
1407 N(2,g) = ux*uy*(ux+uy-3.0)/4.0;
1408 N(3,g) = lx*uy*(lx+uy-3.0)/4.0;
1409 N(4,g) = mx*ly*0.50;
1410 N(5,g) = ux*my*0.50;
1411 N(6,g) = mx*uy*0.50;
1412 N(7,g) = lx*my*0.50;
1413
1414 auto& Nxi = face.Nx;
1415 Nxi(0,0,g) = -ly*(lx+ly-3.0+lx)/4.0;
1416 Nxi(1,0,g) = -lx*(lx+ly-3.0+ly)/4.0;
1417
1418 Nxi(0,1,g) = ly*(ux+ly-3.0+ux)/4.0;
1419 Nxi(1,1,g) = -ux*(ux+ly-3.0+ly)/4.0;
1420
1421 Nxi(0,2,g) = uy*(ux+uy-3.0+ux)/4.0;
1422 Nxi(1,2,g) = ux*(ux+uy-3.0+uy)/4.0;
1423
1424 Nxi(0,3,g) = -uy*(lx+uy-3.0+lx)/4.0;
1425 Nxi(1,3,g) = lx*(lx+uy-3.0+uy)/4.0;
1426
1427 Nxi(0,4,g) = (lx - ux)*ly*0.50;
1428 Nxi(1,4,g) = -mx*0.50;
1429
1430 Nxi(0,5,g) = my*0.50;
1431 Nxi(1,5,g) = (ly - uy)*ux*0.50;
1432
1433 Nxi(0,6,g) = (lx - ux)*uy*0.50;
1434 Nxi(1,6,g) = mx*0.50;
1435
1436 Nxi(0,7,g) = -my*0.50;
1437 Nxi(1,7,g) = (ly - uy)*lx*0.50;
1438 }
1439 },
1440
1441 {ElementType::QUD9, [](int g, faceType& face) -> void
1442 {
1443 auto& xi = face.xi;
1444 double lx = 1.0 - xi(0,g);
1445 double ly = 1.0 - xi(1,g);
1446 double ux = 1.0 + xi(0,g);
1447 double uy = 1.0 + xi(1,g);
1448 double mx = xi(0,g);
1449 double my = xi(1,g);
1450
1451 auto& N = face.N;
1452 N(0,g) = mx*lx*my*ly/4.0;
1453 N(1,g) = -mx*ux*my*ly/4.0;
1454 N(2,g) = mx*ux*my*uy/4.0;
1455 N(3,g) = -mx*lx*my*uy/4.0;
1456 N(4,g) = -lx*ux*my*ly*0.50;
1457 N(5,g) = mx*ux*ly*uy*0.50;
1458 N(6,g) = lx*ux*my*uy*0.50;
1459 N(7,g) = -mx*lx*ly*uy*0.50;
1460 N(8,g) = lx*ux*ly*uy;
1461
1462 auto& Nx = face.Nx;
1463 Nx(0,0,g) = (lx - mx)*my*ly/4.0;
1464 Nx(1,0,g) = (ly - my)*mx*lx/4.0;
1465 Nx(0,1,g) = -(ux + mx)*my*ly/4.0;
1466 Nx(1,1,g) = -(ly - my)*mx*ux/4.0;
1467 Nx(0,2,g) = (ux + mx)*my*uy/4.0;
1468 Nx(1,2,g) = (uy + my)*mx*ux/4.0;
1469 Nx(0,3,g) = -(lx - mx)*my*uy/4.0;
1470 Nx(1,3,g) = -(uy + my)*mx*lx/4.0;
1471 Nx(0,4,g) = -(lx - ux)*my*ly*0.50;
1472 Nx(1,4,g) = -(ly - my)*lx*ux*0.50;
1473 Nx(0,5,g) = (ux + mx)*ly*uy*0.50;
1474 Nx(1,5,g) = (ly - uy)*mx*ux*0.50;
1475 Nx(0,6,g) = (lx - ux)*my*uy*0.50;
1476 Nx(1,6,g) = (uy + my)*lx*ux*0.50;
1477 Nx(0,7,g) = -(lx - mx)*ly*uy*0.50;
1478 Nx(1,7,g) = -(ly - uy)*mx*lx*0.50;
1479 Nx(0,8,g) = (lx - ux)*ly*uy;
1480 Nx(1,8,g) = (ly - uy)*lx*ux;
1481 }
1482 },
1483
1484 {ElementType::LIN1, [](int g, faceType& face) -> void
1485 {
1486 face.N(0,g) = 0.5 * (1.0 - face.xi(0,g));
1487 face.N(1,g) = 0.5 * (1.0 + face.xi(0,g));
1488
1489 face.Nx(0,0,g) = -0.5;
1490 face.Nx(0,1,g) = 0.5;
1491 }
1492 },
1493
1494 {ElementType::LIN2, [](int g, faceType& face) -> void
1495 {
1496 auto& xi = face.xi;
1497 auto& N = face.N;
1498 N(0,g) = -xi(0,g)*(1.0 - xi(0,g))*0.50;
1499 N(1,g) = xi(0,g)*(1.0 + xi(0,g))*0.50;
1500 N(2,g) = (1.0 - xi(0,g))*(1.0 + xi(0,g));
1501
1502 auto& Nx = face.Nx;
1503 Nx(0,0,g) = -0.50 + xi(0,g);
1504 Nx(0,1,g) = 0.50 + xi(0,g);
1505 Nx(0,2,g) = -2.0*xi(0,g);
1506 }
1507 },
1508
1509 {ElementType::QUD4, [](int g, faceType& face) -> void {
1510 auto& xi = face.xi;
1511 double lx = 1.0 - xi(0,g);
1512 double ly = 1.0 - xi(1,g);
1513 double ux = 1.0 + xi(0,g);
1514 double uy = 1.0 + xi(1,g);
1515
1516 auto& N =face.N;
1517 N(0,g) = lx*ly / 4.0;
1518 N(1,g) = ux*ly / 4.0;
1519 N(2,g) = ux*uy / 4.0;
1520 N(3,g) = lx*uy / 4.0;
1521
1522 auto& Nx = face.Nx;
1523 Nx(0,0,g) = -ly / 4.0;
1524 Nx(1,0,g) = -lx / 4.0;
1525 Nx(0,1,g) = ly / 4.0;
1526 Nx(1,1,g) = -ux / 4.0;
1527 Nx(0,2,g) = uy / 4.0;
1528 Nx(1,2,g) = ux / 4.0;
1529 Nx(0,3,g) = -uy / 4.0;
1530 Nx(1,3,g) = lx / 4.0;
1531 }
1532 },
1533
1534 {ElementType::TRI3, [](int g, faceType& face) -> void
1535 {
1536 face.N(0,g) = face.xi(0,g);
1537 face.N(1,g) = face.xi(1,g);
1538 face.N(2,g) = 1.0 - face.xi(0,g) - face.xi(1,g);
1539
1540 face.Nx(0,0,g) = 1.0;
1541 face.Nx(1,0,g) = 0.0;
1542
1543 face.Nx(0,1,g) = 0.0;
1544 face.Nx(1,1,g) = 1.0;
1545
1546 face.Nx(0,2,g) = -1.0;
1547 face.Nx(1,2,g) = -1.0;
1548 }
1549 },
1550
1551 {ElementType::TRI6, [](int g, faceType& face) -> void
1552 {
1553 auto& xi = face.xi;
1554 auto& N = face.N;
1555
1556 double s = 1.0 - xi(0,g) - xi(1,g);
1557 N(0,g) = xi(0,g)*( 2.0*xi(0,g) - 1.0 );
1558 N(1,g) = xi(1,g)*( 2.0*xi(1,g) - 1.0 );
1559 N(2,g) = s *( 2.0*s - 1.0 );
1560 N(3,g) = 4.0*xi(0,g)*xi(1,g);
1561 N(4,g) = 4.0*xi(1,g)*s;
1562 N(5,g) = 4.0*xi(0,g)*s;
1563
1564 auto& Nxi = face.Nx;
1565 Nxi(0,0,g) = 4.0*xi(0,g) - 1.0;
1566 Nxi(1,0,g) = 0.0;
1567
1568 Nxi(0,1,g) = 0.0;
1569 Nxi(1,1,g) = 4.0*xi(1,g) - 1.0;
1570
1571 Nxi(0,2,g) = 1.0 - 4.0*s;
1572 Nxi(1,2,g) = 1.0 - 4.0*s;
1573
1574 Nxi(0,3,g) = 4.0*xi(1,g);
1575 Nxi(1,3,g) = 4.0*xi(0,g);
1576
1577 Nxi(0,4,g) = -4.0*xi(1,g);
1578 Nxi(1,4,g) = 4.0*( s - xi(1,g) );
1579
1580 Nxi(0,5,g) = 4.0*( s - xi(0,g) );
1581 Nxi(1,5,g) = -4.0*xi(0,g);
1582 }
1583 },
1584
1585
1586};
Array3
The Array3 template class implements a simple interface to 3D arrays.
Definition Array3.h:25
faceType
The face type containing mesh at boundary.
Definition ComMod.h:511
mshType
This is the container for a mesh or NURBS patch, those specific to NURBS are noted.
Definition ComMod.h:863
Generated by doxygen 1.9.8