 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
1 /*=========================================================================
2
3 Module: $RCSfile: V_GaussIntegration.cpp,v $
4
5 Copyright (c) 2006 Sandia Corporation.
6 All rights reserved.
7 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
8
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the above copyright notice for more information.
12
13 =========================================================================*/
14
15 /*
16 *
17 * GaussIntegration.cpp performs Gauss Integrations
18 *
19 * This file is part of VERDICT
20 *
21 */
22
23 #define VERDICT_EXPORTS
24
25 #include "moab/verdict.h"
26 #include "V_GaussIntegration.hpp"
27
28 #include <cmath>
29
30 double verdictSqrt2 = sqrt( (double)2.0 );
31
32 int numberGaussPoints;
33 int numberNodes;
34 int numberDims;
35 double gaussPointY[maxNumberGaussPoints];
36 double gaussWeight[maxNumberGaussPoints];
37 double shapeFunction[maxTotalNumberGaussPoints][maxNumberNodes];
38 double dndy1GaussPts[maxTotalNumberGaussPoints][maxNumberNodes];
39 double dndy2GaussPts[maxTotalNumberGaussPoints][maxNumberNodes];
40 double dndy3GaussPts[maxTotalNumberGaussPoints][maxNumberNodes];
41 double totalGaussWeight[maxTotalNumberGaussPoints];
42 int totalNumberGaussPts;
43 double y1Area[maxNumberGaussPointsTri];
44 double y2Area[maxNumberGaussPointsTri];
45 double y1Volume[maxNumberGaussPointsTet];
46 double y2Volume[maxNumberGaussPointsTet];
47 double y3Volume[maxNumberGaussPointsTet];
48 double y4Volume[maxNumberGaussPointsTet];
49
50 void GaussIntegration::initialize( int n, int m, int dim, int tri )
51 {
52 numberGaussPoints = n;
53 numberNodes = m;
54 numberDims = dim;
55
56 if( tri == 1 )
57 // triangular element
58 {
59 if( numberDims == 2 )
60 totalNumberGaussPts = numberGaussPoints;
61 else if( numberDims == 3 )
62 totalNumberGaussPts = numberGaussPoints;
63 }
64 else if( tri == 0 )
65 {
66 if( numberDims == 2 )
67 totalNumberGaussPts = numberGaussPoints * numberGaussPoints;
68 else if( numberDims == 3 )
69 totalNumberGaussPts = numberGaussPoints * numberGaussPoints * numberGaussPoints;
70 }
71 }
72
73 void GaussIntegration::get_shape_func( double shape_function[],
74 double dndy1_at_gauss_pts[],
75 double dndy2_at_gauss_pts[],
76 double gauss_weight[] )
77 {
78 int i, j;
79 for( i = 0; i < totalNumberGaussPts; i++ )
80 {
81 for( j = 0; j < numberNodes; j++ )
82 {
83 shape_function[i * maxNumberNodes + j] = shapeFunction[i][j];
84 dndy1_at_gauss_pts[i * maxNumberNodes + j] = dndy1GaussPts[i][j];
85 dndy2_at_gauss_pts[i * maxNumberNodes + j] = dndy2GaussPts[i][j];
86 }
87 }
88
89 for( i = 0; i < totalNumberGaussPts; i++ )
90 gauss_weight[i] = totalGaussWeight[i];
91 }
92
93 void GaussIntegration::get_shape_func( double shape_function[],
94 double dndy1_at_gauss_pts[],
95 double dndy2_at_gauss_pts[],
96 double dndy3_at_gauss_pts[],
97 double gauss_weight[] )
98 {
99 int i, j;
100 for( i = 0; i < totalNumberGaussPts; i++ )
101 {
102 for( j = 0; j < numberNodes; j++ )
103 {
104 shape_function[i * maxNumberNodes + j] = shapeFunction[i][j];
105 dndy1_at_gauss_pts[i * maxNumberNodes + j] = dndy1GaussPts[i][j];
106 dndy2_at_gauss_pts[i * maxNumberNodes + j] = dndy2GaussPts[i][j];
107 dndy3_at_gauss_pts[i * maxNumberNodes + j] = dndy3GaussPts[i][j];
108 }
109 }
110
111 for( i = 0; i < totalNumberGaussPts; i++ )
112 gauss_weight[i] = totalGaussWeight[i];
113 }
114
115 void GaussIntegration::get_gauss_pts_and_weight()
116 {
117
118 switch( numberGaussPoints )
119 {
120 case 1:
121 gaussPointY[0] = 0.0;
122 gaussWeight[0] = 2.0;
123 break;
124 case 2:
125 gaussPointY[0] = -0.577350269189626;
126 gaussPointY[1] = 0.577350269189626;
127 gaussWeight[0] = 1.0;
128 gaussWeight[1] = 1.0;
129 break;
130 case 3:
131 gaussPointY[0] = -0.774596669241483;
132 gaussPointY[1] = 0.0;
133 gaussPointY[2] = 0.774596669241483;
134 gaussWeight[0] = 0.555555555555555;
135 gaussWeight[1] = 0.888888888888889;
136 gaussWeight[2] = 0.555555555555555;
137 break;
138 }
139 }
140
141 void GaussIntegration::calculate_shape_function_2d_quad()
142 {
143 int ife = 0, i, j;
144 double y1, y2;
145 get_gauss_pts_and_weight();
146
147 switch( numberNodes )
148 {
149 case 4:
150 for( i = 0; i < numberGaussPoints; i++ )
151 {
152 for( j = 0; j < numberGaussPoints; j++ )
153 {
154 y1 = gaussPointY[i];
155 y2 = gaussPointY[j];
156 shapeFunction[ife][0] = 0.25 * ( 1 - y1 ) * ( 1 - y2 );
157 shapeFunction[ife][1] = 0.25 * ( 1 + y1 ) * ( 1 - y2 );
158 shapeFunction[ife][2] = 0.25 * ( 1 + y1 ) * ( 1 + y2 );
159 shapeFunction[ife][3] = 0.25 * ( 1 - y1 ) * ( 1 + y2 );
160
161 dndy1GaussPts[ife][0] = -0.25 * ( 1 - y2 );
162 dndy1GaussPts[ife][1] = 0.25 * ( 1 - y2 );
163 dndy1GaussPts[ife][2] = 0.25 * ( 1 + y2 );
164 dndy1GaussPts[ife][3] = -0.25 * ( 1 + y2 );
165
166 dndy2GaussPts[ife][0] = -0.25 * ( 1 - y1 );
167 dndy2GaussPts[ife][1] = -0.25 * ( 1 + y1 );
168 dndy2GaussPts[ife][2] = 0.25 * ( 1 + y1 );
169 dndy2GaussPts[ife][3] = 0.25 * ( 1 - y1 );
170
171 totalGaussWeight[ife] = gaussWeight[i] * gaussWeight[j];
172 ife++;
173 }
174 }
175 break;
176 case 8:
177 for( i = 0; i < numberGaussPoints; i++ )
178 {
179 for( j = 0; j < numberGaussPoints; j++ )
180 {
181 y1 = gaussPointY[i];
182 y2 = gaussPointY[j];
183 shapeFunction[ife][0] = 0.25 * ( 1 - y1 ) * ( 1 - y2 ) * ( -y1 - y2 - 1 );
184 shapeFunction[ife][1] = 0.25 * ( 1 + y1 ) * ( 1 - y2 ) * ( y1 - y2 - 1 );
185 shapeFunction[ife][2] = 0.25 * ( 1 + y1 ) * ( 1 + y2 ) * ( y1 + y2 - 1 );
186 shapeFunction[ife][3] = 0.25 * ( 1 - y1 ) * ( 1 + y2 ) * ( -y1 + y2 - 1 );
187 shapeFunction[ife][4] = 0.5 * ( 1 - y1 * y1 ) * ( 1 - y2 );
188 shapeFunction[ife][5] = 0.5 * ( 1 - y2 * y2 ) * ( 1 + y1 );
189 shapeFunction[ife][6] = 0.5 * ( 1 - y1 * y1 ) * ( 1 + y2 );
190 shapeFunction[ife][7] = 0.5 * ( 1 - y2 * y2 ) * ( 1 - y1 );
191
192 dndy1GaussPts[ife][0] = 0.25 * ( 1 - y2 ) * ( 2.0 * y1 + y2 );
193 dndy1GaussPts[ife][1] = 0.25 * ( 1 - y2 ) * ( 2.0 * y1 - y2 );
194 dndy1GaussPts[ife][2] = 0.25 * ( 1 + y2 ) * ( 2.0 * y1 + y2 );
195 dndy1GaussPts[ife][3] = 0.25 * ( 1 + y2 ) * ( 2.0 * y1 - y2 );
196
197 dndy1GaussPts[ife][4] = -y1 * ( 1 - y2 );
198 dndy1GaussPts[ife][5] = 0.5 * ( 1 - y2 * y2 );
199 dndy1GaussPts[ife][6] = -y1 * ( 1 + y2 );
200 dndy1GaussPts[ife][7] = -0.5 * ( 1 - y2 * y2 );
201
202 dndy2GaussPts[ife][0] = 0.25 * ( 1 - y1 ) * ( 2.0 * y2 + y1 );
203 dndy2GaussPts[ife][1] = 0.25 * ( 1 + y1 ) * ( 2.0 * y2 - y1 );
204 dndy2GaussPts[ife][2] = 0.25 * ( 1 + y1 ) * ( 2.0 * y2 + y1 );
205 dndy2GaussPts[ife][3] = 0.25 * ( 1 - y1 ) * ( 2.0 * y2 - y1 );
206
207 dndy2GaussPts[ife][4] = -0.5 * ( 1 - y1 * y1 );
208 dndy2GaussPts[ife][5] = -y2 * ( 1 + y1 );
209 dndy2GaussPts[ife][6] = 0.5 * ( 1 - y1 * y1 );
210 dndy2GaussPts[ife][7] = -y2 * ( 1 - y1 );
211
212 totalGaussWeight[ife] = gaussWeight[i] * gaussWeight[j];
213 ife++;
214 }
215 }
216 break;
217 }
218 }
219
220 void GaussIntegration::calculate_shape_function_3d_hex()
221 {
222 int ife = 0, i, j, k, node_id;
223 double y1, y2, y3, sign_node_y1, sign_node_y2, sign_node_y3;
224 double y1_term, y2_term, y3_term, y123_temp;
225
226 get_gauss_pts_and_weight();
227
228 switch( numberNodes )
229 {
230 case 8:
231 for( i = 0; i < numberGaussPoints; i++ )
232 {
233 for( j = 0; j < numberGaussPoints; j++ )
234 {
235 for( k = 0; k < numberGaussPoints; k++ )
236 {
237 y1 = gaussPointY[i];
238 y2 = gaussPointY[j];
239 y3 = gaussPointY[k];
240
241 for( node_id = 0; node_id < numberNodes; node_id++ )
242 {
243 get_signs_for_node_local_coord_hex( node_id, sign_node_y1, sign_node_y2, sign_node_y3 );
244
245 y1_term = 1 + sign_node_y1 * y1;
246 y2_term = 1 + sign_node_y2 * y2;
247 y3_term = 1 + sign_node_y3 * y3;
248
249 shapeFunction[ife][node_id] = 0.125 * y1_term * y2_term * y3_term;
250 dndy1GaussPts[ife][node_id] = 0.125 * sign_node_y1 * y2_term * y3_term;
251 dndy2GaussPts[ife][node_id] = 0.125 * sign_node_y2 * y1_term * y3_term;
252 dndy3GaussPts[ife][node_id] = 0.125 * sign_node_y3 * y1_term * y2_term;
253 }
254 totalGaussWeight[ife] = gaussWeight[i] * gaussWeight[j] * gaussWeight[k];
255 ife++;
256 }
257 }
258 }
259 break;
260 case 20:
261 for( i = 0; i < numberGaussPoints; i++ )
262 {
263 for( j = 0; j < numberGaussPoints; j++ )
264 {
265 for( k = 0; k < numberGaussPoints; k++ )
266 {
267 y1 = gaussPointY[i];
268 y2 = gaussPointY[j];
269 y3 = gaussPointY[k];
270
271 for( node_id = 0; node_id < numberNodes; node_id++ )
272 {
273 get_signs_for_node_local_coord_hex( node_id, sign_node_y1, sign_node_y2, sign_node_y3 );
274
275 y1_term = 1 + sign_node_y1 * y1;
276 y2_term = 1 + sign_node_y2 * y2;
277 y3_term = 1 + sign_node_y3 * y3;
278 y123_temp = sign_node_y1 * y1 + sign_node_y2 * y2 + sign_node_y3 * y3 - 2.;
279
280 switch( node_id )
281 {
282 case 0:
283 case 1:
284 case 2:
285 case 3:
286 case 4:
287 case 5:
288 case 6:
289 case 7: {
290 shapeFunction[ife][node_id] = 0.125 * y1_term * y2_term * y3_term * y123_temp;
291 dndy1GaussPts[ife][node_id] = 0.125 * sign_node_y1 * y123_temp * y2_term * y3_term +
292 0.125 * y1_term * y2_term * y3_term * sign_node_y1;
293 dndy2GaussPts[ife][node_id] = 0.125 * sign_node_y2 * y1_term * y3_term * y123_temp +
294 0.125 * y1_term * y2_term * y3_term * sign_node_y2;
295 dndy3GaussPts[ife][node_id] = 0.125 * sign_node_y3 * y1_term * y2_term * y123_temp +
296 0.125 * y1_term * y2_term * y3_term * sign_node_y3;
297 break;
298 }
299 case 8:
300 case 10:
301 case 16:
302 case 18: {
303 shapeFunction[ife][node_id] = 0.25 * ( 1 - y1 * y1 ) * y2_term * y3_term;
304 dndy1GaussPts[ife][node_id] = -0.5 * y1 * y2_term * y3_term;
305 dndy2GaussPts[ife][node_id] = 0.25 * ( 1 - y1 * y1 ) * sign_node_y2 * y3_term;
306 dndy3GaussPts[ife][node_id] = 0.25 * ( 1 - y1 * y1 ) * y2_term * sign_node_y3;
307 break;
308 }
309 case 9:
310 case 11:
311 case 17:
312 case 19: {
313 shapeFunction[ife][node_id] = 0.25 * ( 1 - y2 * y2 ) * y1_term * y3_term;
314 dndy1GaussPts[ife][node_id] = 0.25 * ( 1 - y2 * y2 ) * sign_node_y1 * y3_term;
315 dndy2GaussPts[ife][node_id] = -0.5 * y2 * y1_term * y3_term;
316 dndy3GaussPts[ife][node_id] = 0.25 * ( 1 - y2 * y2 ) * y1_term * sign_node_y3;
317 break;
318 }
319 case 12:
320 case 13:
321 case 14:
322 case 15: {
323 shapeFunction[ife][node_id] = 0.25 * ( 1 - y3 * y3 ) * y1_term * y2_term;
324 dndy1GaussPts[ife][node_id] = 0.25 * ( 1 - y3 * y3 ) * sign_node_y1 * y2_term;
325 dndy2GaussPts[ife][node_id] = 0.25 * ( 1 - y3 * y3 ) * y1_term * sign_node_y2;
326 dndy3GaussPts[ife][node_id] = -0.5 * y3 * y1_term * y2_term;
327 break;
328 }
329 }
330 }
331 totalGaussWeight[ife] = gaussWeight[i] * gaussWeight[j] * gaussWeight[k];
332 ife++;
333 }
334 }
335 }
336 break;
337 }
338 }
339
340 void GaussIntegration::calculate_derivative_at_nodes( double dndy1_at_nodes[][maxNumberNodes],
341 double dndy2_at_nodes[][maxNumberNodes] )
342 {
343 double y1 = 0., y2 = 0.;
344 int i;
345 for( i = 0; i < numberNodes; i++ )
346 {
347 switch( i )
348 {
349 case 0:
350 y1 = -1.;
351 y2 = -1.;
352 break;
353 case 1:
354 y1 = 1.;
355 y2 = -1.;
356 break;
357 case 2:
358 y1 = 1.;
359 y2 = 1.;
360 break;
361 case 3:
362 y1 = -1.;
363 y2 = 1.;
364 break;
365
366 // midside nodes if there is any
367
368 case 4:
369 y1 = 0.;
370 y2 = -1.;
371 break;
372
373 case 5:
374 y1 = 1.;
375 y2 = 0.;
376 break;
377
378 case 6:
379 y1 = 0.;
380 y2 = 1.;
381 break;
382
383 case 7:
384 y1 = -1.;
385 y2 = 0.;
386 break;
387 }
388
389 switch( numberNodes )
390 {
391 case 4:
392 // dn_i/dy1 evaluated at node i
393 dndy1_at_nodes[i][0] = -0.25 * ( 1 - y2 );
394 dndy1_at_nodes[i][1] = 0.25 * ( 1 - y2 );
395 dndy1_at_nodes[i][2] = 0.25 * ( 1 + y2 );
396 dndy1_at_nodes[i][3] = -0.25 * ( 1 + y2 );
397
398 // dn_i/dy2 evaluated at node i
399 dndy2_at_nodes[i][0] = -0.25 * ( 1 - y1 );
400 dndy2_at_nodes[i][1] = -0.25 * ( 1 + y1 );
401 dndy2_at_nodes[i][2] = 0.25 * ( 1 + y1 );
402 dndy2_at_nodes[i][3] = 0.25 * ( 1 - y1 );
403 break;
404
405 case 8:
406
407 dndy1_at_nodes[i][0] = 0.25 * ( 1 - y2 ) * ( 2.0 * y1 + y2 );
408 dndy1_at_nodes[i][1] = 0.25 * ( 1 - y2 ) * ( 2.0 * y1 - y2 );
409 dndy1_at_nodes[i][2] = 0.25 * ( 1 + y2 ) * ( 2.0 * y1 + y2 );
410 dndy1_at_nodes[i][3] = 0.25 * ( 1 + y2 ) * ( 2.0 * y1 - y2 );
411
412 dndy1_at_nodes[i][4] = -y1 * ( 1 - y2 );
413 dndy1_at_nodes[i][5] = 0.5 * ( 1 - y2 * y2 );
414 dndy1_at_nodes[i][6] = -y1 * ( 1 + y2 );
415 dndy1_at_nodes[i][7] = -0.5 * ( 1 - y2 * y2 );
416
417 dndy2_at_nodes[i][0] = 0.25 * ( 1 - y1 ) * ( 2.0 * y2 + y1 );
418 dndy2_at_nodes[i][1] = 0.25 * ( 1 + y1 ) * ( 2.0 * y2 - y1 );
419 dndy2_at_nodes[i][2] = 0.25 * ( 1 + y1 ) * ( 2.0 * y2 + y1 );
420 dndy2_at_nodes[i][3] = 0.25 * ( 1 - y1 ) * ( 2.0 * y2 - y1 );
421
422 dndy2_at_nodes[i][4] = -0.5 * ( 1 - y1 * y1 );
423 dndy2_at_nodes[i][5] = -y2 * ( 1 + y1 );
424 dndy2_at_nodes[i][6] = 0.5 * ( 1 - y1 * y1 );
425 dndy2_at_nodes[i][7] = -y2 * ( 1 - y1 );
426 break;
427 }
428 }
429 }
430
431 void GaussIntegration::calculate_derivative_at_nodes_3d( double dndy1_at_nodes[][maxNumberNodes],
432 double dndy2_at_nodes[][maxNumberNodes],
433 double dndy3_at_nodes[][maxNumberNodes] )
434 {
435 double y1, y2, y3, sign_node_y1, sign_node_y2, sign_node_y3;
436 double y1_term, y2_term, y3_term, y123_temp;
437 int node_id, node_id_2;
438 for( node_id = 0; node_id < numberNodes; node_id++ )
439 {
440 get_signs_for_node_local_coord_hex( node_id, y1, y2, y3 );
441
442 switch( numberNodes )
443 {
444 case 8:
445 for( node_id_2 = 0; node_id_2 < numberNodes; node_id_2++ )
446 {
447 get_signs_for_node_local_coord_hex( node_id_2, sign_node_y1, sign_node_y2, sign_node_y3 );
448 y1_term = 1 + sign_node_y1 * y1;
449 y2_term = 1 + sign_node_y2 * y2;
450 y3_term = 1 + sign_node_y3 * y3;
451
452 dndy1_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y1 * y2_term * y3_term;
453
454 dndy2_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y2 * y1_term * y3_term;
455
456 dndy3_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y3 * y1_term * y2_term;
457 }
458 break;
459 case 20:
460 for( node_id_2 = 0; node_id_2 < numberNodes; node_id_2++ )
461 {
462 get_signs_for_node_local_coord_hex( node_id_2, sign_node_y1, sign_node_y2, sign_node_y3 );
463
464 y1_term = 1 + sign_node_y1 * y1;
465 y2_term = 1 + sign_node_y2 * y2;
466 y3_term = 1 + sign_node_y3 * y3;
467 y123_temp = sign_node_y1 * y1 + sign_node_y2 * y2 + sign_node_y3 * y3 - 2.;
468 switch( node_id_2 )
469 {
470 case 0:
471 case 1:
472 case 2:
473 case 3:
474 case 4:
475 case 5:
476 case 6:
477 case 7: {
478 dndy1_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y1 * y2_term * y3_term * y123_temp +
479 0.125 * y1_term * y2_term * y3_term * sign_node_y1;
480 dndy2_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y2 * y1_term * y3_term * y123_temp +
481 0.125 * y1_term * y2_term * y3_term * sign_node_y2;
482 dndy3_at_nodes[node_id][node_id_2] = 0.125 * sign_node_y3 * y1_term * y2_term * y123_temp +
483 0.125 * y1_term * y2_term * y3_term * sign_node_y3;
484 break;
485 }
486 case 8:
487 case 10:
488 case 16:
489 case 18: {
490 dndy1_at_nodes[node_id][node_id_2] = -0.5 * y1 * y2_term * y3_term;
491 dndy2_at_nodes[node_id][node_id_2] = 0.25 * ( 1 - y1 * y1 ) * sign_node_y2 * y3_term;
492 dndy3_at_nodes[node_id][node_id_2] = 0.25 * ( 1 - y1 * y1 ) * y2_term * sign_node_y3;
493 break;
494 }
495 case 9:
496 case 11:
497 case 17:
498 case 19: {
499 dndy1_at_nodes[node_id][node_id_2] = 0.25 * ( 1 - y2 * y2 ) * sign_node_y1 * y3_term;
500 dndy2_at_nodes[node_id][node_id_2] = -0.5 * y2 * y1_term * y3_term;
501 dndy3_at_nodes[node_id][node_id_2] = 0.25 * ( 1 - y2 * y2 ) * y1_term * sign_node_y3;
502 break;
503 }
504 case 12:
505 case 13:
506 case 14:
507 case 15: {
508 dndy1_at_nodes[node_id][node_id_2] = 0.25 * ( 1 - y3 * y3 ) * sign_node_y1 * y2_term;
509 dndy2_at_nodes[node_id][node_id_2] = 0.25 * ( 1 - y3 * y3 ) * y1_term * sign_node_y2;
510 dndy3_at_nodes[node_id][node_id_2] = -0.5 * y3 * y1_term * y2_term;
511 break;
512 }
513 }
514 }
515 break;
516 }
517 }
518 }
519
520 void GaussIntegration::get_signs_for_node_local_coord_hex( int node_id,
521 double& sign_node_y1,
522 double& sign_node_y2,
523 double& sign_node_y3 )
524 {
525 switch( node_id )
526 {
527 case 0:
528 sign_node_y1 = -1.;
529 sign_node_y2 = -1.;
530 sign_node_y3 = -1.;
531 break;
532 case 1:
533 sign_node_y1 = 1.;
534 sign_node_y2 = -1.;
535 sign_node_y3 = -1.;
536 break;
537 case 2:
538 sign_node_y1 = 1.;
539 sign_node_y2 = 1.;
540 sign_node_y3 = -1.;
541 break;
542 case 3:
543 sign_node_y1 = -1.;
544 sign_node_y2 = 1.;
545 sign_node_y3 = -1.;
546 break;
547 case 4:
548 sign_node_y1 = -1.;
549 sign_node_y2 = -1.;
550 sign_node_y3 = 1.;
551 break;
552 case 5:
553 sign_node_y1 = 1.;
554 sign_node_y2 = -1.;
555 sign_node_y3 = 1.;
556 break;
557 case 6:
558 sign_node_y1 = 1.;
559 sign_node_y2 = 1.;
560 sign_node_y3 = 1.;
561 break;
562 case 7:
563 sign_node_y1 = -1.;
564 sign_node_y2 = 1.;
565 sign_node_y3 = 1.;
566 break;
567 case 8:
568 sign_node_y1 = 0;
569 sign_node_y2 = -1.;
570 sign_node_y3 = -1.;
571 break;
572 case 9:
573 sign_node_y1 = 1.;
574 sign_node_y2 = 0;
575 sign_node_y3 = -1.;
576 break;
577 case 10:
578 sign_node_y1 = 0;
579 sign_node_y2 = 1.;
580 sign_node_y3 = -1.;
581 break;
582 case 11:
583 sign_node_y1 = -1.;
584 sign_node_y2 = 0.;
585 sign_node_y3 = -1.;
586 break;
587 case 12:
588 sign_node_y1 = -1.;
589 sign_node_y2 = -1.;
590 sign_node_y3 = 0.;
591 break;
592 case 13:
593 sign_node_y1 = 1.;
594 sign_node_y2 = -1.;
595 sign_node_y3 = 0.;
596 break;
597 case 14:
598 sign_node_y1 = 1.;
599 sign_node_y2 = 1.;
600 sign_node_y3 = 0.;
601 break;
602 case 15:
603 sign_node_y1 = -1.;
604 sign_node_y2 = 1.;
605 sign_node_y3 = 0.;
606 break;
607 case 16:
608 sign_node_y1 = 0;
609 sign_node_y2 = -1.;
610 sign_node_y3 = 1.;
611 break;
612 case 17:
613 sign_node_y1 = 1.;
614 sign_node_y2 = 0;
615 sign_node_y3 = 1.;
616 break;
617 case 18:
618 sign_node_y1 = 0;
619 sign_node_y2 = 1.;
620 sign_node_y3 = 1.;
621 break;
622 case 19:
623 sign_node_y1 = -1.;
624 sign_node_y2 = 0.;
625 sign_node_y3 = 1.;
626 break;
627 }
628 }
629
630 void GaussIntegration::get_tri_rule_pts_and_weight()
631 {
632 // get triangular rule integration points and weight
633
634 switch( numberGaussPoints )
635 {
636 case 6:
637 y1Area[0] = 0.09157621;
638 y2Area[0] = 0.09157621;
639
640 y1Area[1] = 0.09157621;
641 y2Area[1] = 0.8168476;
642
643 y1Area[2] = 0.8168476;
644 y2Area[2] = 0.09157621;
645
646 y1Area[3] = 0.4459485;
647 y2Area[3] = 0.4459485;
648
649 y1Area[4] = 0.4459485;
650 y2Area[4] = 0.1081030;
651
652 y1Area[5] = 0.1081030;
653 y2Area[5] = 0.4459485;
654
655 int i;
656 for( i = 0; i < 3; i++ )
657 {
658 totalGaussWeight[i] = 0.06348067;
659 }
660 for( i = 3; i < 6; i++ )
661 {
662 totalGaussWeight[i] = 0.1289694;
663 }
664 break;
665 }
666 }
667
668 void GaussIntegration::calculate_shape_function_2d_tri()
669 {
670 int ife;
671 double y1, y2, y3;
672 get_tri_rule_pts_and_weight();
673
674 for( ife = 0; ife < totalNumberGaussPts; ife++ )
675 {
676 y1 = y1Area[ife];
677 y2 = y2Area[ife];
678 y3 = 1.0 - y1 - y2;
679
680 shapeFunction[ife][0] = y1 * ( 2. * y1 - 1. );
681 shapeFunction[ife][1] = y2 * ( 2. * y2 - 1. );
682 shapeFunction[ife][2] = y3 * ( 2. * y3 - 1. );
683
684 shapeFunction[ife][3] = 4. * y1 * y2;
685 shapeFunction[ife][4] = 4. * y2 * y3;
686 shapeFunction[ife][5] = 4. * y1 * y3;
687
688 dndy1GaussPts[ife][0] = 4 * y1 - 1.;
689 dndy1GaussPts[ife][1] = 0;
690 dndy1GaussPts[ife][2] = 1 - 4. * y3;
691
692 dndy1GaussPts[ife][3] = 4. * y2;
693 dndy1GaussPts[ife][4] = -4. * y2;
694 dndy1GaussPts[ife][5] = 4. * ( 1 - 2 * y1 - y2 );
695
696 dndy2GaussPts[ife][0] = 0.0;
697 dndy2GaussPts[ife][1] = 4. * y2 - 1.;
698 dndy2GaussPts[ife][2] = 1 - 4. * y3;
699
700 dndy2GaussPts[ife][3] = 4. * y1;
701 dndy2GaussPts[ife][4] = 4. * ( 1 - y1 - 2. * y2 );
702 dndy2GaussPts[ife][5] = -4. * y1;
703 }
704 }
705
706 void GaussIntegration::calculate_derivative_at_nodes_2d_tri( double dndy1_at_nodes[][maxNumberNodes],
707 double dndy2_at_nodes[][maxNumberNodes] )
708 {
709 double y1 = 0., y2 = 0., y3;
710 int i;
711 for( i = 0; i < numberNodes; i++ )
712 {
713 switch( i )
714 {
715 case 0:
716 y1 = 1.;
717 y2 = 0.;
718 break;
719 case 1:
720 y1 = 0.;
721 y2 = 1.;
722 break;
723 case 2:
724 y1 = 0.;
725 y2 = 0.;
726 break;
727 case 3:
728 y1 = 0.5;
729 y2 = 0.5;
730 break;
731 case 4:
732 y1 = 0.;
733 y2 = 0.5;
734 break;
735 case 5:
736 y1 = 0.5;
737 y2 = 0.0;
738 break;
739 }
740
741 y3 = 1. - y1 - y2;
742
743 dndy1_at_nodes[i][0] = 4 * y1 - 1.;
744 dndy1_at_nodes[i][1] = 0;
745 dndy1_at_nodes[i][2] = 1 - 4. * y3;
746
747 dndy1_at_nodes[i][3] = 4. * y2;
748 dndy1_at_nodes[i][4] = -4. * y2;
749 dndy1_at_nodes[i][5] = 4. * ( 1 - 2 * y1 - y2 );
750
751 dndy2_at_nodes[i][0] = 0.0;
752 dndy2_at_nodes[i][1] = 4. * y2 - 1.;
753 dndy2_at_nodes[i][2] = 1 - 4. * y3;
754
755 dndy2_at_nodes[i][3] = 4. * y1;
756 dndy2_at_nodes[i][4] = 4. * ( 1 - y1 - 2. * y2 );
757 dndy2_at_nodes[i][5] = -4. * y1;
758 }
759 }
760 void GaussIntegration::get_tet_rule_pts_and_weight()
761 {
762 // get tetrahedron rule integration points and weight
763
764 double a, b;
765 switch( numberGaussPoints )
766 {
767 case 1:
768 // 1 integration point formula, degree of precision 1
769 y1Volume[0] = 0.25;
770 y2Volume[0] = 0.25;
771 y3Volume[0] = 0.25;
772 y4Volume[0] = 0.25;
773 totalGaussWeight[0] = 1.;
774 break;
775 case 4:
776 // 4 integration points formula, degree of precision 2
777 a = 0.58541020;
778 b = 0.13819660;
779
780 y1Volume[0] = a;
781 y2Volume[0] = b;
782 y3Volume[0] = b;
783 y4Volume[0] = b;
784
785 y1Volume[1] = b;
786 y2Volume[1] = a;
787 y3Volume[1] = b;
788 y4Volume[1] = b;
789
790 y1Volume[2] = b;
791 y2Volume[2] = b;
792 y3Volume[2] = a;
793 y4Volume[2] = b;
794
795 y1Volume[3] = b;
796 y2Volume[3] = b;
797 y3Volume[3] = b;
798 y4Volume[3] = a;
799
800 int i;
801 for( i = 0; i < 4; i++ )
802 {
803 totalGaussWeight[i] = 0.25;
804 }
805 break;
806 }
807 }
808
809 void GaussIntegration::calculate_shape_function_3d_tet()
810 {
811 int ife;
812 double y1, y2, y3, y4;
813 get_tet_rule_pts_and_weight();
814
815 switch( numberNodes )
816 {
817 case 10: // 10 nodes quadratic tet
818 {
819 for( ife = 0; ife < totalNumberGaussPts; ife++ )
820 {
821 // y1,y2,y3,y4 are the volume coordinates
822 y1 = y1Volume[ife];
823 y2 = y2Volume[ife];
824 y3 = y3Volume[ife];
825 y4 = y4Volume[ife];
826
827 // shape function is the same as in ABAQUS
828 // it is different from that in all the FEA book
829 // in which node is the first node
830 // here at node 1 y4=1
831 shapeFunction[ife][0] = y4 * ( 2. * y4 - 1. );
832 shapeFunction[ife][1] = y1 * ( 2. * y1 - 1. );
833 shapeFunction[ife][2] = y2 * ( 2. * y2 - 1. );
834 shapeFunction[ife][3] = y3 * ( 2. * y3 - 1. );
835
836 shapeFunction[ife][4] = 4. * y1 * y4;
837 shapeFunction[ife][5] = 4. * y1 * y2;
838 shapeFunction[ife][6] = 4. * y2 * y4;
839 shapeFunction[ife][7] = 4. * y3 * y4;
840 shapeFunction[ife][8] = 4. * y1 * y3;
841 shapeFunction[ife][9] = 4. * y2 * y3;
842
843 dndy1GaussPts[ife][0] = 1 - 4 * y4;
844 dndy1GaussPts[ife][1] = 4 * y1 - 1.;
845 dndy1GaussPts[ife][2] = 0;
846 dndy1GaussPts[ife][3] = 0;
847
848 dndy1GaussPts[ife][4] = 4. * ( y4 - y1 );
849 dndy1GaussPts[ife][5] = 4. * y2;
850 dndy1GaussPts[ife][6] = -4. * y2;
851 dndy1GaussPts[ife][7] = -4. * y3;
852 dndy1GaussPts[ife][8] = 4. * y3;
853 dndy1GaussPts[ife][9] = 0;
854
855 dndy2GaussPts[ife][0] = 1 - 4 * y4;
856 dndy2GaussPts[ife][1] = 0;
857 dndy2GaussPts[ife][2] = 4. * y2 - 1.;
858 dndy2GaussPts[ife][3] = 0;
859
860 dndy2GaussPts[ife][4] = -4. * y1;
861 dndy2GaussPts[ife][5] = 4. * y1;
862 dndy2GaussPts[ife][6] = 4. * ( y4 - y2 );
863 dndy2GaussPts[ife][7] = -4. * y3;
864 dndy2GaussPts[ife][8] = 0.;
865 dndy2GaussPts[ife][9] = 4. * y3;
866
867 dndy3GaussPts[ife][0] = 1 - 4 * y4;
868 dndy3GaussPts[ife][1] = 0;
869 dndy3GaussPts[ife][2] = 0;
870 dndy3GaussPts[ife][3] = 4. * y3 - 1.;
871
872 dndy3GaussPts[ife][4] = -4. * y1;
873 dndy3GaussPts[ife][5] = 0;
874 dndy3GaussPts[ife][6] = -4. * y2;
875 dndy3GaussPts[ife][7] = 4. * ( y4 - y3 );
876 dndy3GaussPts[ife][8] = 4. * y1;
877 dndy3GaussPts[ife][9] = 4. * y2;
878 }
879 break;
880 }
881 case 4: // four node linear tet for debug purpose
882 {
883 for( ife = 0; ife < totalNumberGaussPts; ife++ )
884 {
885 y1 = y1Volume[ife];
886 y2 = y2Volume[ife];
887 y3 = y3Volume[ife];
888 y4 = y4Volume[ife];
889
890 shapeFunction[ife][0] = y4;
891 shapeFunction[ife][1] = y1;
892 shapeFunction[ife][2] = y2;
893 shapeFunction[ife][3] = y3;
894
895 dndy1GaussPts[ife][0] = -1.;
896 dndy1GaussPts[ife][1] = 1;
897 dndy1GaussPts[ife][2] = 0;
898 dndy1GaussPts[ife][3] = 0;
899
900 dndy2GaussPts[ife][0] = -1.;
901 dndy2GaussPts[ife][1] = 0;
902 dndy2GaussPts[ife][2] = 1;
903 dndy2GaussPts[ife][3] = 0;
904
905 dndy3GaussPts[ife][0] = -1.;
906 dndy3GaussPts[ife][1] = 0;
907 dndy3GaussPts[ife][2] = 0;
908 dndy3GaussPts[ife][3] = 1;
909 }
910 break;
911 }
912 }
913 }
914
915 void GaussIntegration::calculate_derivative_at_nodes_3d_tet( double dndy1_at_nodes[][maxNumberNodes],
916 double dndy2_at_nodes[][maxNumberNodes],
917 double dndy3_at_nodes[][maxNumberNodes] )
918 {
919 double y1, y2, y3, y4;
920 int i;
921
922 switch( numberNodes )
923 {
924 case 10: {
925 for( i = 0; i < numberNodes; i++ )
926 {
927 get_node_local_coord_tet( i, y1, y2, y3, y4 );
928
929 dndy1_at_nodes[i][0] = 1 - 4 * y4;
930 dndy1_at_nodes[i][1] = 4 * y1 - 1.;
931 dndy1_at_nodes[i][2] = 0;
932 dndy1_at_nodes[i][3] = 0;
933
934 dndy1_at_nodes[i][4] = 4. * ( y4 - y1 );
935 dndy1_at_nodes[i][5] = 4. * y2;
936 dndy1_at_nodes[i][6] = -4. * y2;
937 dndy1_at_nodes[i][7] = -4. * y3;
938 dndy1_at_nodes[i][8] = 4. * y3;
939 dndy1_at_nodes[i][9] = 0;
940
941 dndy2_at_nodes[i][0] = 1 - 4 * y4;
942 dndy2_at_nodes[i][1] = 0;
943 dndy2_at_nodes[i][2] = 4. * y2 - 1.;
944 dndy2_at_nodes[i][3] = 0;
945 dndy2_at_nodes[i][4] = -4. * y1;
946 dndy2_at_nodes[i][5] = 4. * y1;
947 dndy2_at_nodes[i][6] = 4. * ( y4 - y2 );
948 dndy2_at_nodes[i][7] = -4. * y3;
949 dndy2_at_nodes[i][8] = 0.;
950 dndy2_at_nodes[i][9] = 4. * y3;
951
952 dndy3_at_nodes[i][0] = 1 - 4 * y4;
953 dndy3_at_nodes[i][1] = 0;
954 dndy3_at_nodes[i][2] = 0;
955 dndy3_at_nodes[i][3] = 4. * y3 - 1.;
956
957 dndy3_at_nodes[i][4] = -4. * y1;
958 dndy3_at_nodes[i][5] = 0;
959 dndy3_at_nodes[i][6] = -4. * y2;
960 dndy3_at_nodes[i][7] = 4. * ( y4 - y3 );
961 dndy3_at_nodes[i][8] = 4. * y1;
962 dndy3_at_nodes[i][9] = 4. * y2;
963 }
964 break;
965 }
966 case 4: {
967 for( i = 0; i < numberNodes; i++ )
968 {
969 get_node_local_coord_tet( i, y1, y2, y3, y4 );
970 dndy1_at_nodes[i][0] = -1.;
971 dndy1_at_nodes[i][1] = 1;
972 dndy1_at_nodes[i][2] = 0;
973 dndy1_at_nodes[i][3] = 0;
974
975 dndy2_at_nodes[i][0] = -1.;
976 dndy2_at_nodes[i][1] = 0;
977 dndy2_at_nodes[i][2] = 1;
978 dndy2_at_nodes[i][3] = 0;
979
980 dndy3_at_nodes[i][0] = -1.;
981 dndy3_at_nodes[i][1] = 0;
982 dndy3_at_nodes[i][2] = 0;
983 dndy3_at_nodes[i][3] = 1;
984 }
985 break;
986 }
987 }
988 }
989
990 void GaussIntegration::get_node_local_coord_tet( int node_id, double& y1, double& y2, double& y3, double& y4 )
991 {
992 switch( node_id )
993 {
994 case 0:
995 y1 = 0.;
996 y2 = 0.;
997 y3 = 0.;
998 y4 = 1.;
999 break;
1000 case 1:
1001 y1 = 1.;
1002 y2 = 0.;
1003 y3 = 0.;
1004 y4 = 0.;
1005 break;
1006 case 2:
1007 y1 = 0.;
1008 y2 = 1.;
1009 y3 = 0.;
1010 y4 = 0.;
1011 break;
1012 case 3:
1013 y1 = 0.;
1014 y2 = 0.;
1015 y3 = 1.;
1016 y4 = 0.;
1017 break;
1018 case 4:
1019 y1 = 0.5;
1020 y2 = 0.;
1021 y3 = 0.;
1022 y4 = 0.5;
1023 break;
1024 case 5:
1025 y1 = 0.5;
1026 y2 = 0.5;
1027 y3 = 0.;
1028 y4 = 0.;
1029 break;
1030 case 6:
1031 y1 = 0.;
1032 y2 = 0.5;
1033 y3 = 0.;
1034 y4 = 0.5;
1035 break;
1036 case 7:
1037 y1 = 0.;
1038 y2 = 0.0;
1039 y3 = 0.5;
1040 y4 = 0.5;
1041 break;
1042 case 8:
1043 y1 = 0.5;
1044 y2 = 0.;
1045 y3 = 0.5;
1046 y4 = 0.0;
1047 break;
1048 case 9:
1049 y1 = 0.;
1050 y2 = 0.5;
1051 y3 = 0.5;
1052 y4 = 0.;
1053 break;
1054 }
1055 }
1.9.1.