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