svMultiPhysics
Loading...
Searching...
No Matches
nn_elem_gip.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 Gauss integration data.
5//
6using GetElementGausIntMapType = std::map<ElementType, std::function<void(const int, const int,
7 Vector<double>&, Array<double>&)>>;
8
9/// @brief Replicates 'SUBROUTINE GETGIP(insd, eType, nG, w, xi)' defined in NN.f.
10//
11GetElementGausIntMapType get_element_gauss_int_data = {
12
13 {ElementType::HEX8, [](const int insd, const int nG, Vector<double>& w, Array<double>& xi) -> void {
14 w = 1.0;
15 double s = 1.0 / sqrt(3.0);
16 double t = -1.0 / sqrt(3.0);
17
18 xi(0,0) = t;
19 xi(1,0) = t;
20 xi(2,0) = t;
21
22 xi(0,1) = s;
23 xi(1,1) = t;
24 xi(2,1) = t;
25
26 xi(0,2) = s;
27 xi(1,2) = s;
28 xi(2,2) = t;
29
30 xi(0,3) = t;
31 xi(1,3) = s;
32 xi(2,3) = t;
33
34 xi(0,4) = t;
35 xi(1,4) = t;
36 xi(2,4) = s;
37
38 xi(0,5) = s;
39 xi(1,5) = t;
40 xi(2,5) = s;
41
42 xi(0,6) = s;
43 xi(1,6) = s;
44 xi(2,6) = s;
45
46 xi(0,7) = t;
47 xi(1,7) = s;
48 xi(2,7) = s;
49 }
50 },
51
52 {ElementType::HEX20, [](const int insd, const int nG, Vector<double>& w, Array<double>& xi) -> void {
53 for (int i = 0; i < 8; i++) {
54 w(i) = 125.0 / 729.0;
55 }
56
57 for (int i = 8; i < 20; i++) {
58 w(i) = 200.0 / 729.0;
59 }
60
61 for (int i = 20; i < 26; i++) {
62 w(i) = 320.0 / 729.0;
63 }
64
65 w(26) = 512.0 / 729.0;
66
67 double s = sqrt(0.6);
68 double t = 0.0;
69
70 int n = 0;
71 xi(0, n) = -s; xi(1, n) = -s; xi(2, n) = -s; n++;
72 xi(0, n) = s; xi(1, n) = -s; xi(2, n) = -s; n++;
73 xi(0, n) = s; xi(1, n) = s; xi(2, n) = -s; n++;
74 xi(0, n) = -s; xi(1, n) = s; xi(2, n) = -s; n++;
75 xi(0, n) = -s; xi(1, n) = -s; xi(2, n) = s; n++;
76 xi(0, n) = s; xi(1, n) = -s; xi(2, n) = s; n++;
77 xi(0, n) = s; xi(1, n) = s; xi(2, n) = s; n++;
78 xi(0, n) = -s; xi(1, n) = s; xi(2, n) = s; n++;
79
80 xi(0, n) = t; xi(1, n) = -s; xi(2, n) = -s; n++;
81 xi(0,n) = s; xi(1,n) = t; xi(2,n) = -s; n++;
82 xi(0,n) = t; xi(1,n) = s; xi(2,n) = -s; n++;
83 xi(0,n) = -s; xi(1,n) = t; xi(2,n) = -s; n++;
84 xi(0,n) = t; xi(1,n) = -s; xi(2,n) = s; n++;
85 xi(0,n) = s; xi(1,n) = t; xi(2,n) = s; n++;
86 xi(0,n) = t; xi(1,n) = s; xi(2,n) = s; n++;
87 xi(0,n) = -s; xi(1,n) = t; xi(2,n) = s; n++;
88 xi(0,n) = -s; xi(1,n) = -s; xi(2,n) = t; n++;
89 xi(0,n) = s; xi(1,n) = -s; xi(2,n) = t; n++;
90 xi(0,n) = s; xi(1,n) = s; xi(2,n) = t; n++;
91 xi(0,n) = -s; xi(1,n) = s; xi(2,n) = t; n++;
92
93 xi(0,n) = -s; xi(1,n) = t; xi(2,n) = t; n++;
94 xi(0,n) = s; xi(1,n) = t; xi(2,n) = t; n++;
95 xi(0,n) = t; xi(1,n) = -s; xi(2,n) = t; n++;
96 xi(0,n) = t; xi(1,n) = s; xi(2,n) = t; n++;
97 xi(0,n) = t; xi(1,n) = t; xi(2,n) = -s; n++;
98 xi(0,n) = t; xi(1,n) = t; xi(2,n) = s; n++;
99
100 xi(0,n) = t; xi(1,n) = t; xi(2,n) = t;
101 }
102 },
103
104 {ElementType::HEX27, [](const int insd, const int nG, Vector<double>& w, Array<double>& xi) -> void {
105
106 for (int i = 0; i < 8; i++) {
107 w(i) = 125.0 / 729.0;
108 }
109
110 for (int i = 8; i < 20; i++) {
111 w(i) = 200.0 / 729.0;
112 }
113
114 for (int i = 20; i < 26; i++) {
115 w(i) = 320.0 / 729.0;
116 }
117
118 w(26) = 512.0 / 729.0;
119
120 double s = sqrt(0.6);
121 double t = 0.0;
122
123 int n = 0;
124 xi(0, n) = -s; xi(1, n) = -s; xi(2, n) = -s; n++;
125 xi(0, n) = s; xi(1, n) = -s; xi(2, n) = -s; n++;
126 xi(0, n) = s; xi(1, n) = s; xi(2, n) = -s; n++;
127 xi(0, n) = -s; xi(1, n) = s; xi(2, n) = -s; n++;
128 xi(0, n) = -s; xi(1, n) = -s; xi(2, n) = s; n++;
129 xi(0, n) = s; xi(1, n) = -s; xi(2, n) = s; n++;
130 xi(0, n) = s; xi(1, n) = s; xi(2, n) = s; n++;
131 xi(0, n) = -s; xi(1, n) = s; xi(2, n) = s; n++;
132
133 xi(0, n) = t; xi(1, n) = -s; xi(2, n) = -s; n++;
134 xi(0,n) = s; xi(1,n) = t; xi(2,n) = -s; n++;
135 xi(0,n) = t; xi(1,n) = s; xi(2,n) = -s; n++;
136 xi(0,n) = -s; xi(1,n) = t; xi(2,n) = -s; n++;
137 xi(0,n) = t; xi(1,n) = -s; xi(2,n) = s; n++;
138 xi(0,n) = s; xi(1,n) = t; xi(2,n) = s; n++;
139 xi(0,n) = t; xi(1,n) = s; xi(2,n) = s; n++;
140 xi(0,n) = -s; xi(1,n) = t; xi(2,n) = s; n++;
141 xi(0,n) = -s; xi(1,n) = -s; xi(2,n) = t; n++;
142 xi(0,n) = s; xi(1,n) = -s; xi(2,n) = t; n++;
143 xi(0,n) = s; xi(1,n) = s; xi(2,n) = t; n++;
144 xi(0,n) = -s; xi(1,n) = s; xi(2,n) = t; n++;
145
146 xi(0,n) = -s; xi(1,n) = t; xi(2,n) = t; n++;
147 xi(0,n) = s; xi(1,n) = t; xi(2,n) = t; n++;
148 xi(0,n) = t; xi(1,n) = -s; xi(2,n) = t; n++;
149 xi(0,n) = t; xi(1,n) = s; xi(2,n) = t; n++;
150 xi(0,n) = t; xi(1,n) = t; xi(2,n) = -s; n++;
151 xi(0,n) = t; xi(1,n) = t; xi(2,n) = s; n++;
152
153 xi(0,n) = t; xi(1,n) = t; xi(2,n) = t;
154 }
155 },
156
157 {ElementType::LIN1, [](const int insd, const int nG, Vector<double>& w, Array<double>& xi) -> void {
158 w = 1.0;
159 double s = 1.0 / sqrt(3.0);
160 xi(0,0) = -s;
161 xi(0,1) = s;
162 }
163 },
164
165 {ElementType::QUD4, [](const int insd, const int nG, Vector<double>& w, Array<double>& xi) -> void {
166 w = 1.0;
167 double s = 1.0 / sqrt(3.0);
168 xi(0,0) = -s; xi(1,0) = -s;
169 xi(0,1) = s; xi(1,1) = -s;
170 xi(0,2) = s; xi(1,2) = s;
171 xi(0,3) = -s; xi(1,3) = s;
172 }
173 },
174
175 {ElementType::QUD9, [](const int insd, const int nG, Vector<double>& w, Array<double>& xi) -> void {
176 for (int i = 0; i < 4; i++) {
177 w(i) = 25.0 / 81.0;
178 w(i+4) = 40.0 / 81.0;
179 }
180 w(8) = 64.0 / 81.0;
181
182 double s = sqrt(6.0);
183 double t = 0.0;
184
185 xi(0,0) = -s;
186 xi(1,0) = -s;
187
188 xi(0,1) = s;
189 xi(1,1) = -s;
190
191 xi(0,2) = s;
192 xi(1,2) = s;
193
194 xi(0,3) = -s;
195 xi(1,3) = s;
196
197 xi(0,4) = t;
198 xi(1,4) = -s;
199
200 xi(0,5) = s;
201 xi(1,5) = t;
202
203 xi(0,6) = t;
204 xi(1,6) = s;
205
206 xi(0,7) = -s;
207 xi(1,7) = t;
208
209 xi(0,8) = t;
210 xi(1,8) = t;
211 }
212 },
213
214 {ElementType::TET4, [](const int insd, const int nG, Vector<double>& w, Array<double>& xi) -> void {
215 for (int i = 0; i < nG; i++) {
216 w(i) = 1.0 / 24.0;
217 }
218 double s = (5.0 + 3.0*sqrt(5.0)) / 20.0;
219 double t = (5.0 - sqrt(5.0)) / 20.0;
220 xi(0,0) = s; xi(1,0) = t; xi(2,0) = t;
221 xi(0,1) = t; xi(1,1) = s; xi(2,1) = t;
222 xi(0,2) = t; xi(1,2) = t; xi(2,2) = s;
223 xi(0,3) = t; xi(1,3) = t; xi(2,3) = t;
224 }
225 },
226
227 {ElementType::TRI3, [](const int insd, const int nG, Vector<double>& w, Array<double>& xi) -> void {
228 for (int i = 0; i < nG; i++) {
229 w(i) = 1.0 / 6.0;
230 }
231 double s = 2.0/3.0;
232 double t = 1.0/6.0;
233 xi(0,0) = t;
234 xi(1,0) = t;
235 xi(0,1) = s;
236 xi(1,1) = t;
237 xi(0,2) = t;
238 xi(1,2) = s;
239 }
240 },
241
242 {ElementType::WDG, [](const int insd, const int nG, Vector<double>& w, Array<double>& xi) -> void {
243 for (int i = 0; i < nG; i++) {
244 w(i) = 1.0 / 6.0;
245 }
246 double s = 2.0 / 3.0;
247 double t = 1.0 / 6.0;
248 double uz = 1.0 / sqrt(3.0);
249 double lz = -1.0 / sqrt(3.0);
250 xi(0,0) = s; xi(1,0) = t; xi(2,0) = lz;
251 xi(0,1) = t; xi(1,1) = s; xi(2,1) = lz;
252 xi(0,2) = t; xi(1,2) = t; xi(2,2) = lz;
253 xi(0,3) = s; xi(1,3) = t; xi(2,3) = uz;
254 xi(0,4) = t; xi(1,4) = s; xi(2,4) = uz;
255 xi(0,5) = t; xi(1,5) = t; xi(2,5) = uz;
256 }
257 },
258
259};
260
261
262/// @brief Define a map type used to set element Gauss integration data.
263using SetElementGausIntMapType = std::map<ElementType, std::function<void(mshType&)>>;
264
265/// @brief Replicates 'SUBROUTINE GETGIP(insd, eType, nG, w, xi)' defined in NN.f.
266///
267/// \todo [TODO:DaveP] I screwed up here, should just have a single map for mesh and face types.
268//
269SetElementGausIntMapType set_element_gauss_int_data = {
270
271 {ElementType::HEX8, [](mshType& mesh) -> void {
272 for (int i = 0; i < mesh.nG; i++) {
273 mesh.w(i) = 1.0;
274 }
275 double s = 1.0 / sqrt(3.0);
276 double t = -1.0 / sqrt(3.0);
277 mesh.xi(0,0) = t; mesh.xi(1,0) = t; mesh.xi(2,0) = t;
278 mesh.xi(0,1) = s; mesh.xi(1,1) = t; mesh.xi(2,1) = t;
279 mesh.xi(0,2) = s; mesh.xi(1,2) = s; mesh.xi(2,2) = t;
280 mesh.xi(0,3) = t; mesh.xi(1,3) = s; mesh.xi(2,3) = t;
281 mesh.xi(0,4) = t; mesh.xi(1,4) = t; mesh.xi(2,4) = s;
282 mesh.xi(0,5) = s; mesh.xi(1,5) = t; mesh.xi(2,5) = s;
283 mesh.xi(0,6) = s; mesh.xi(1,6) = s; mesh.xi(2,6) = s;
284 mesh.xi(0,7) = t; mesh.xi(1,7) = s; mesh.xi(2,7) = s;
285 }
286 },
287
288 {ElementType::HEX20, [](mshType& mesh) -> void {
289
290 for (int i = 0; i < 8; i++) {
291 mesh.w(i) = 125.0 / 729.0;
292 }
293
294 for (int i = 8; i < 20; i++) {
295 mesh.w(i) = 200.0 / 729.0;
296 }
297
298 for (int i = 20; i < 26; i++) {
299 mesh.w(i) = 320.0 / 729.0;
300 }
301
302 mesh.w(26) = 512.0 / 729.0;
303
304 double s = sqrt(0.6);
305 double t = 0.0;
306
307 int n = 0;
308 mesh.xi(0, n) = -s; mesh.xi(1, n) = -s; mesh.xi(2, n) = -s; n++;
309 mesh.xi(0, n) = s; mesh.xi(1, n) = -s; mesh.xi(2, n) = -s; n++;
310 mesh.xi(0, n) = s; mesh.xi(1, n) = s; mesh.xi(2, n) = -s; n++;
311 mesh.xi(0, n) = -s; mesh.xi(1, n) = s; mesh.xi(2, n) = -s; n++;
312 mesh.xi(0, n) = -s; mesh.xi(1, n) = -s; mesh.xi(2, n) = s; n++;
313 mesh.xi(0, n) = s; mesh.xi(1, n) = -s; mesh.xi(2, n) = s; n++;
314 mesh.xi(0, n) = s; mesh.xi(1, n) = s; mesh.xi(2, n) = s; n++;
315 mesh.xi(0, n) = -s; mesh.xi(1, n) = s; mesh.xi(2, n) = s; n++;
316
317 mesh.xi(0, n) = t; mesh.xi(1, n) = -s; mesh.xi(2, n) = -s; n++;
318 mesh.xi(0,n) = s; mesh.xi(1,n) = t; mesh.xi(2,n) = -s; n++;
319 mesh.xi(0,n) = t; mesh.xi(1,n) = s; mesh.xi(2,n) = -s; n++;
320 mesh.xi(0,n) = -s; mesh.xi(1,n) = t; mesh.xi(2,n) = -s; n++;
321 mesh.xi(0,n) = t; mesh.xi(1,n) = -s; mesh.xi(2,n) = s; n++;
322 mesh.xi(0,n) = s; mesh.xi(1,n) = t; mesh.xi(2,n) = s; n++;
323 mesh.xi(0,n) = t; mesh.xi(1,n) = s; mesh.xi(2,n) = s; n++;
324 mesh.xi(0,n) = -s; mesh.xi(1,n) = t; mesh.xi(2,n) = s; n++;
325 mesh.xi(0,n) = -s; mesh.xi(1,n) = -s; mesh.xi(2,n) = t; n++;
326 mesh.xi(0,n) = s; mesh.xi(1,n) = -s; mesh.xi(2,n) = t; n++;
327 mesh.xi(0,n) = s; mesh.xi(1,n) = s; mesh.xi(2,n) = t; n++;
328 mesh.xi(0,n) = -s; mesh.xi(1,n) = s; mesh.xi(2,n) = t; n++;
329
330 mesh.xi(0,n) = -s; mesh.xi(1,n) = t; mesh.xi(2,n) = t; n++;
331 mesh.xi(0,n) = s; mesh.xi(1,n) = t; mesh.xi(2,n) = t; n++;
332 mesh.xi(0,n) = t; mesh.xi(1,n) = -s; mesh.xi(2,n) = t; n++;
333 mesh.xi(0,n) = t; mesh.xi(1,n) = s; mesh.xi(2,n) = t; n++;
334 mesh.xi(0,n) = t; mesh.xi(1,n) = t; mesh.xi(2,n) = -s; n++;
335 mesh.xi(0,n) = t; mesh.xi(1,n) = t; mesh.xi(2,n) = s; n++;
336
337 mesh.xi(0,n) = t; mesh.xi(1,n) = t; mesh.xi(2,n) = t;
338 }
339 },
340
341 {ElementType::HEX27, [](mshType& mesh) -> void {
342
343 for (int i = 0; i < 8; i++) {
344 mesh.w(i) = 125.0 / 729.0;
345 }
346
347 for (int i = 8; i < 20; i++) {
348 mesh.w(i) = 200.0 / 729.0;
349 }
350
351 for (int i = 20; i < 26; i++) {
352 mesh.w(i) = 320.0 / 729.0;
353 }
354
355 mesh.w(26) = 512.0 / 729.0;
356
357 double s = sqrt(0.6);
358 double t = 0.0;
359
360 int n = 0;
361 mesh.xi(0, n) = -s; mesh.xi(1, n) = -s; mesh.xi(2, n) = -s; n++;
362 mesh.xi(0, n) = s; mesh.xi(1, n) = -s; mesh.xi(2, n) = -s; n++;
363 mesh.xi(0, n) = s; mesh.xi(1, n) = s; mesh.xi(2, n) = -s; n++;
364 mesh.xi(0, n) = -s; mesh.xi(1, n) = s; mesh.xi(2, n) = -s; n++;
365 mesh.xi(0, n) = -s; mesh.xi(1, n) = -s; mesh.xi(2, n) = s; n++;
366 mesh.xi(0, n) = s; mesh.xi(1, n) = -s; mesh.xi(2, n) = s; n++;
367 mesh.xi(0, n) = s; mesh.xi(1, n) = s; mesh.xi(2, n) = s; n++;
368 mesh.xi(0, n) = -s; mesh.xi(1, n) = s; mesh.xi(2, n) = s; n++;
369
370 mesh.xi(0, n) = t; mesh.xi(1, n) = -s; mesh.xi(2, n) = -s; n++;
371 mesh.xi(0,n) = s; mesh.xi(1,n) = t; mesh.xi(2,n) = -s; n++;
372 mesh.xi(0,n) = t; mesh.xi(1,n) = s; mesh.xi(2,n) = -s; n++;
373 mesh.xi(0,n) = -s; mesh.xi(1,n) = t; mesh.xi(2,n) = -s; n++;
374 mesh.xi(0,n) = t; mesh.xi(1,n) = -s; mesh.xi(2,n) = s; n++;
375 mesh.xi(0,n) = s; mesh.xi(1,n) = t; mesh.xi(2,n) = s; n++;
376 mesh.xi(0,n) = t; mesh.xi(1,n) = s; mesh.xi(2,n) = s; n++;
377 mesh.xi(0,n) = -s; mesh.xi(1,n) = t; mesh.xi(2,n) = s; n++;
378 mesh.xi(0,n) = -s; mesh.xi(1,n) = -s; mesh.xi(2,n) = t; n++;
379 mesh.xi(0,n) = s; mesh.xi(1,n) = -s; mesh.xi(2,n) = t; n++;
380 mesh.xi(0,n) = s; mesh.xi(1,n) = s; mesh.xi(2,n) = t; n++;
381 mesh.xi(0,n) = -s; mesh.xi(1,n) = s; mesh.xi(2,n) = t; n++;
382
383 mesh.xi(0,n) = -s; mesh.xi(1,n) = t; mesh.xi(2,n) = t; n++;
384 mesh.xi(0,n) = s; mesh.xi(1,n) = t; mesh.xi(2,n) = t; n++;
385 mesh.xi(0,n) = t; mesh.xi(1,n) = -s; mesh.xi(2,n) = t; n++;
386 mesh.xi(0,n) = t; mesh.xi(1,n) = s; mesh.xi(2,n) = t; n++;
387 mesh.xi(0,n) = t; mesh.xi(1,n) = t; mesh.xi(2,n) = -s; n++;
388 mesh.xi(0,n) = t; mesh.xi(1,n) = t; mesh.xi(2,n) = s; n++;
389
390 mesh.xi(0,n) = t; mesh.xi(1,n) = t; mesh.xi(2,n) = t;
391 }
392 },
393
394 {ElementType::LIN1, [](mshType& mesh) -> void {
395 for (int i = 0; i < mesh.nG; i++) {
396 mesh.w(i) = 1.0;
397 }
398 double s = 1.0 / sqrt(3.0);
399 mesh.xi(0,0) = -s;
400 mesh.xi(0,1) = s;
401 }
402 },
403
404 {ElementType::LIN2, [](mshType& mesh) -> void {
405 mesh.w[0] = 5.0 / 9.0;
406 mesh.w[1] = 5.0 / 9.0;
407 mesh.w[2] = 8.0 / 9.0;
408
409 double s = sqrt(0.6);
410
411 mesh.xi(0,0) = -s;
412 mesh.xi(0,1) = s;
413 mesh.xi(0,2) = 0.0;
414 }
415 },
416
417 {ElementType::QUD4, [](mshType& mesh) -> void {
418 for (int i = 0; i < mesh.nG; i++) {
419 mesh.w[i] = 1.0;
420 }
421 double s = 1.0 / sqrt(3.0);
422 mesh.xi(0,0) = -s; mesh.xi(1,0) = -s;
423 mesh.xi(0,1) = s; mesh.xi(1,1) = -s;
424 mesh.xi(0,2) = s; mesh.xi(1,2) = s;
425 mesh.xi(0,3) = -s; mesh.xi(1,3) = s;
426 }
427 },
428
429 {ElementType::QUD9, [](mshType& mesh) -> void {
430 mesh.w(0) = 25.0 / 81.0;
431 mesh.w(1) = 25.0 / 81.0;
432 mesh.w(2) = 25.0 / 81.0;
433 mesh.w(3) = 25.0 / 81.0;
434
435 mesh.w(4) = 40.0 / 81.0;
436 mesh.w(5) = 40.0 / 81.0;
437 mesh.w(6) = 40.0 / 81.0;
438 mesh.w(7) = 40.0 / 81.0;
439 mesh.w(8) = 64.0 / 81.0;
440
441 double s = sqrt(0.6);
442 double t = 0.0;
443
444 mesh.xi(0,0) = -s;
445 mesh.xi(1,0) = -s;
446
447 mesh.xi(0,1) = s;
448 mesh.xi(1,1) = -s;
449
450 mesh.xi(0,2) = s;
451 mesh.xi(1,2) = s;
452
453 mesh.xi(0,3) = -s;
454 mesh.xi(1,3) = s;
455
456 mesh.xi(0,4) = t;
457 mesh.xi(1,4) = -s;
458
459 mesh.xi(0,5) = s;
460 mesh.xi(1,5) = t;
461
462 mesh.xi(0,6) = t;
463 mesh.xi(1,6) = s;
464
465 mesh.xi(0,7) = -s;
466 mesh.xi(1,7) = t;
467
468 mesh.xi(0,8) = t;
469 mesh.xi(1,8) = t;
470 }
471 },
472
473
474 {ElementType::TET4, [](mshType& mesh) -> void {
475 for (int i = 0; i < mesh.nG; i++) {
476 mesh.w(i) = 1.0 / 24.0;
477 }
478 // s=0.25: central quadrature rule
479 // s=(5.0+3.0*sqrt(5.0))/20.0: original gaussian quadrature rule
480 // s=1.0: nodal quadrature
481 double s = mesh.qmTET4;
482 // t=0.25: central quadrature rule
483 // t=(5.0-sqrt(5.0))/20.0: original gaussian quadrature rule
484 // t=0.0: nodal quadrature
485 double t = (1-s)/3.0;
486 mesh.xi(0,0) = s; mesh.xi(1,0) = t; mesh.xi(2,0) = t;
487 mesh.xi(0,1) = t; mesh.xi(1,1) = s; mesh.xi(2,1) = t;
488 mesh.xi(0,2) = t; mesh.xi(1,2) = t; mesh.xi(2,2) = s;
489 mesh.xi(0,3) = t; mesh.xi(1,3) = t; mesh.xi(2,3) = t;
490 }
491 },
492
493 {ElementType::TET10, [](mshType& mesh) -> void {
494 mesh.w(0) = 0.0302836780970890;
495
496 mesh.w(1) = 0.0060267857142860;
497 mesh.w(2) = 0.0060267857142860;
498 mesh.w(3) = 0.0060267857142860;
499 mesh.w(4) = 0.0060267857142860;
500
501 mesh.w(5) = 0.0116452490860290;
502 mesh.w(6) = 0.0116452490860290;
503 mesh.w(7) = 0.0116452490860290;
504 mesh.w(8) = 0.0116452490860290;
505
506 mesh.w(9) = 0.0109491415613860;
507 mesh.w(10) = 0.0109491415613860;
508 mesh.w(11) = 0.0109491415613860;
509 mesh.w(12) = 0.0109491415613860;
510 mesh.w(13) = 0.0109491415613860;
511 mesh.w(14) = 0.0109491415613860;
512
513 double s = 0.250;
514 auto& xi = mesh.xi;
515 xi(0,0) = s; xi(1,0) = s; xi(2,0) = s;
516
517 s = 0.3333333333333330;
518 double t = 0.0;
519 xi(0,1) = t; xi(1,1) = s; xi(2,1) = s;
520 xi(0,2) = s; xi(1,2) = t; xi(2,2) = s;
521 xi(0,3) = s; xi(1,3) = s; xi(2,3) = t;
522 xi(0,4) = s; xi(1,4) = s; xi(2,4) = s;
523
524 s = 0.0909090909090910;
525 t = 0.7272727272727270;
526 xi(0,5) = t; xi(1,5) = s; xi(2,5) = s;
527 xi(0,6) = s; xi(1,6) = t; xi(2,6) = s;
528 xi(0,7) = s; xi(1,7) = s; xi(2,7) = t;
529 xi(0,8) = s; xi(1,8) = s; xi(2,8) = s;
530
531 s = 0.0665501535736640;
532 t = 0.4334498464263360;
533 xi(0, 9) = s; xi(1, 9) = s; xi(2, 9) = t;
534 xi(0,10) = s; xi(1,10) = t; xi(2,10) = s;
535 xi(0,11) = s; xi(1,11) = t; xi(2,11) = t;
536 xi(0,12) = t; xi(1,12) = t; xi(2,12) = s;
537 xi(0,13) = t; xi(1,13) = s; xi(2,13) = t;
538 xi(0,14) = t; xi(1,14) = s; xi(2,14) = s;
539 }
540 },
541
542 {ElementType::TRI3, [](mshType& mesh) -> void {
543 for (int i = 0; i < mesh.nG; i++) {
544 mesh.w(i) = 1.0 / 6.0;
545 }
546 double s = 2.0 / 3.0;
547 double t = 1.0 / 6.0;
548 mesh.xi(0,0) = t; mesh.xi(1,0) = t;
549 mesh.xi(0,1) = s; mesh.xi(1,1) = t;
550 mesh.xi(0,2) = t; mesh.xi(1,2) = s;
551 }
552 },
553
554 {ElementType::TRI6, [](mshType& mesh) -> void {
555 mesh.w[0] = 0.225000000000000 * 5.0e-1 ;
556 mesh.w[1] = 0.125939180544827 * 5.0e-1 ;
557 mesh.w[2] = 0.125939180544827 * 5.0e-1 ;
558 mesh.w[3] = 0.125939180544827 * 5.0e-1 ;
559 mesh.w(4) = 0.132394152788506 * 5.0e-1 ;
560 mesh.w(5) = 0.132394152788506 * 5.0e-1 ;
561 mesh.w(6) = 0.132394152788506 * 5.0e-1 ;
562
563 double s = 0.333333333333333;
564 mesh.xi(0,0) = s; mesh.xi(1,0) = s;
565
566 s = 0.797426985353087;
567 double t = 0.101286507323456;
568 mesh.xi(0,1) = s; mesh.xi(1,1) = t;
569 mesh.xi(0,2) = t; mesh.xi(1,2) = s;
570 mesh.xi(0,3) = t; mesh.xi(1,3) = t;
571
572 s = 0.059715871789770;
573 t = 0.470142064105115;
574 mesh.xi(0,4) = s; mesh.xi(1,4) = t;
575 mesh.xi(0,5) = t; mesh.xi(1,5) = s;
576 mesh.xi(0,6) = t; mesh.xi(1,6) = t;
577 }
578 },
579
580 {ElementType::WDG, [](mshType& mesh) -> void {
581 for (int i = 0; i < mesh.nG; i++) {
582 mesh.w(i) = 1.0 / 6.0;
583 }
584 double s = 2.0 / 3.0;
585 double t = 1.0 / 6.0;
586 double uz = 1.0 / sqrt(3.0);
587 double lz = -1.0 / sqrt(3.0);
588 mesh.xi(0,0) = s; mesh.xi(1,0) = t; mesh.xi(2,0) = lz;
589 mesh.xi(0,1) = t; mesh.xi(1,1) = s; mesh.xi(2,1) = lz;
590 mesh.xi(0,2) = t; mesh.xi(1,2) = t; mesh.xi(2,2) = lz;
591 mesh.xi(0,3) = s; mesh.xi(1,3) = t; mesh.xi(2,3) = uz;
592 mesh.xi(0,4) = t; mesh.xi(1,4) = s; mesh.xi(2,4) = uz;
593 mesh.xi(0,5) = t; mesh.xi(1,5) = t; mesh.xi(2,5) = uz;
594 }
595 },
596
597};
598
599/// @brief Define a map type used to set face Gauss integration data.
600//
601using SetFaceGausIntMapType = std::map<ElementType, std::function<void(faceType&)>>;
602
603/// \todo [TODO:DaveP] I screwed up here, should just have a single map for mesh and face types.
604//
605SetFaceGausIntMapType set_face_gauss_int_data = {
606
607 {ElementType::PNT, [](faceType& face) -> void {
608 for (int i = 0; i < face.nG; i++) {
609 face.w(i) = 1.0;
610 }
611 }
612 },
613
614 {ElementType::LIN1, [](faceType& face) -> void {
615 for (int i = 0; i < face.nG; i++) {
616 face.w(i) = 1.0;
617 }
618 double s = 1.0 / sqrt(3.0);
619 face.xi(0,0) = -s;
620 face.xi(0,1) = s;
621 }
622 },
623
624 {ElementType::LIN2, [](faceType& face) -> void {
625 face.w[0] = 5.0 / 9.0;
626 face.w[1] = 5.0 / 9.0;
627 face.w[2] = 8.0 / 9.0;
628
629 double s = sqrt(0.6);
630
631 face.xi(0,0) = -s;
632 face.xi(0,1) = s;
633 face.xi(0,2) = 0.0;
634 }
635 },
636
637 {ElementType::QUD4, [](faceType& face) -> void {
638 face.w = 1.0;
639 double s = 1.0 / sqrt(3.0);
640 face.xi(0,0) = -s; face.xi(1,0) = -s;
641 face.xi(0,1) = s; face.xi(1,1) = -s;
642 face.xi(0,2) = s; face.xi(1,2) = s;
643 face.xi(0,3) = -s; face.xi(1,3) = s;
644 }
645 },
646
647 {ElementType::QUD8, [](faceType& face) -> void {
648 face.w(0) = 25.0/81.0;
649 face.w(1) = 25.0/81.0;
650 face.w(2) = 25.0/81.0;
651 face.w(3) = 25.0/81.0;
652
653 face.w(4) = 40.0/81.0;
654 face.w(5) = 40.0/81.0;
655 face.w(6) = 40.0/81.0;
656 face.w(7) = 40.0/81.0;
657
658 face.w(8) = 64.0/81.0;
659
660 double s = sqrt(0.6);
661 double t = 0.0;
662
663 face.xi(0,0) = -s; face.xi(1,0) = -s;
664 face.xi(0,1) = s; face.xi(1,1) = -s;
665 face.xi(0,2) = s; face.xi(1,2) = s;
666 face.xi(0,3) = -s; face.xi(1,3) = s;
667 face.xi(0,4) = t; face.xi(1,4) = -s;
668 face.xi(0,5) = s; face.xi(1,5) = t;
669 face.xi(0,6) = t; face.xi(1,6) = s;
670 face.xi(0,7) = -s; face.xi(1,7) = t;
671 face.xi(0,8) = t; face.xi(1,8) = t;
672 }
673 },
674
675 {ElementType::QUD9, [](faceType& face) -> void {
676 face.w(0) = 25.0 / 81.0;
677 face.w(1) = 25.0 / 81.0;
678 face.w(2) = 25.0 / 81.0;
679 face.w(3) = 25.0 / 81.0;
680
681 face.w(4) = 40.0 / 81.0;
682 face.w(5) = 40.0 / 81.0;
683 face.w(6) = 40.0 / 81.0;
684 face.w(7) = 40.0 / 81.0;
685 face.w(8) = 64.0 / 81.0;
686
687 double s = sqrt(0.6);
688 double t = 0.0;
689
690 face.xi(0,0) = -s;
691 face.xi(1,0) = -s;
692
693 face.xi(0,1) = s;
694 face.xi(1,1) = -s;
695
696 face.xi(0,2) = s;
697 face.xi(1,2) = s;
698
699 face.xi(0,3) = -s;
700 face.xi(1,3) = s;
701
702 face.xi(0,4) = t;
703 face.xi(1,4) = -s;
704
705 face.xi(0,5) = s;
706 face.xi(1,5) = t;
707
708 face.xi(0,6) = t;
709 face.xi(1,6) = s;
710
711 face.xi(0,7) = -s;
712 face.xi(1,7) = t;
713
714 face.xi(0,8) = t;
715 face.xi(1,8) = t;
716 }
717 },
718
719
720 {ElementType::TRI3, [](faceType& face) -> void {
721 //for (int i = 0; i < face.nG; i++) {
722 // face.w(i) = 1.0 / 6.0;
723 //}
724 face.w = 1.0 / 6.0;
725
726 // s=1.0/3.0: central quadrature rule
727 // s=2.0/3.0: original gaussian quadrature rule
728 // s=1.0: nodal quadrature
729 double s = face.qmTRI3;
730 // t=1.0/3.0: central quadrature rule
731 // t=1.0/6.0: original gaussian quadrature rule
732 // t=0.0: nodal quadrature
733 double t = -0.5*s+0.5;
734 face.xi(0,0) = t; face.xi(1,0) = t;
735 face.xi(0,1) = s; face.xi(1,1) = t;
736 face.xi(0,2) = t; face.xi(1,2) = s;
737 }
738 },
739
740 {ElementType::TRI6, [](faceType& face) -> void {
741 face.w[0] = 0.225000000000000 * 5.0e-1;
742 face.w[1] = 0.125939180544827 * 5.0e-1;
743 face.w[2] = 0.125939180544827 * 5.0e-1;
744 face.w[3] = 0.125939180544827 * 5.0e-1 ;
745 face.w(4) = 0.132394152788506 * 5.0e-1 ;
746 face.w(5) = 0.132394152788506 * 5.0e-1 ;
747 face.w(6) = 0.132394152788506 * 5.0e-1 ;
748
749 double s = 0.333333333333333;
750 face.xi(0,0) = s; face.xi(1,0) = s;
751
752 s = 0.797426985353087;
753 double t = 0.101286507323456;
754 face.xi(0,1) = s; face.xi(1,1) = t;
755 face.xi(0,2) = t; face.xi(1,2) = s;
756 face.xi(0,3) = t; face.xi(1,3) = t;
757
758 s = 0.059715871789770;
759 t = 0.470142064105115;
760 face.xi(0,4) = s; face.xi(1,4) = t;
761 face.xi(0,5) = t; face.xi(1,5) = s;
762 face.xi(0,6) = t; face.xi(1,6) = t;
763 }
764 },
765
766
767
768
769};
Vector
The Vector template class is used for storing int and double data.
Definition Vector.h:23
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