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