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