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