Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
V_GaussIntegration.cpp
Go to the documentation of this file.
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 }