FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
static_LIB_shell.f90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
6 use hecmw, only : kint, kreal
8
9 !--------------------------------------------------------------------
10
11 implicit none
12
13 !--------------------------------------------------------------------
14
15 !--------------------------------------------------------------
16 !
17 ! (Programmer)
18 ! Gaku Hashimoto
19 ! Department of Human and Engineered Environmental Studies
20 ! Graduate School of Frontier Sciences, The University of Tokyo
21 ! 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8563 JAPAN
22 !
23 ! (Ref.)
24 ! [1] Noguchi, H. and Hisada, T.,
25 ! "Sensitivity analysis in post-buckling problems of shell
26 ! structures,"
27 ! Computers & Structures, Vol.47, No.4, pp.699-710, (1993).
28 ! [2] Dvorkin, E.N. and Bathe, K.J.,
29 ! "A Continuum Mechanics Based Four-node Shell Element for
30 ! General Non-linear Analysis,"
31 ! Engineering Computations, Vol.1, pp.77-88, (1984).
32 ! [3] Bucalem, M.L. and Bathe, K.J.,
33 ! "Higher-order MITC general shell element,"
34 ! International Journal for Numerical Methods in
35 ! Engineering, Vol.36, pp.3729-3754, (1993).
36 ! [4] Lee, P.S. and Bathe, K.J.,
37 ! "Development of MITC Isotropic Triangular Shell Finite
38 ! Elements,"
39 ! Computers & Structures, Vol.82, pp.945-962, (2004).
40 !
41 ! Xi YUAN
42 ! Apr. 13, 2019: Introduce mass matrix calculation
43 ! (Ref.)
44 ! [5] E.Hinton, T.A.Rock, O.C.Zienkiewicz(1976): A Note on Mass Lumping
45 ! and Related Process in FEM. International Journal on Earthquake Eng
46 ! and structural dynamics, 4, pp245-249
47 !
48 !--------------------------------------------------------------
49
50 !--------------------------------------------------------------------
51
52contains
53
54
55 !####################################################################
56 subroutine stf_shell_mitc &
57 (etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag, nddisp)
58 !####################################################################
59
60 use mmechgauss
62 use m_matmatrix
63
64 !--------------------------------------------------------------------
65
66 integer(kind = kint), intent(in) :: etype
67 integer(kind = kint), intent(in) :: nn, mixflag
68 integer(kind = kint), intent(in) :: ndof
69 real(kind = kreal), intent(in) :: ecoord(3, nn)
70 type(tgaussstatus), intent(in) :: gausses(:)
71 real(kind = kreal), intent(out) :: stiff(:, :)
72 real(kind = kreal), intent(in) :: thick
73
74 real(kind = kreal), intent(in), optional :: nddisp(3, nn)
75
76 !--------------------------------------------------------------------
77
78 integer :: flag, flag_dof
79 integer :: i, j, m
80 integer :: lx, ly
81 integer :: fetype
82 integer :: ny
83 integer :: ntying
84 integer :: npoints_tying(3)
85 integer :: it, ip
86 integer :: na, nb
87 integer :: isize, jsize
88 integer :: jsize1, jsize2, jsize3, &
89 jsize4, jsize5, jsize6
90 integer :: n_layer,n_totlyr, sstable(24)
91
92 real(kind = kreal) :: d(5, 5), b(5, ndof*nn), db(5, ndof*nn)
93 real(kind = kreal) :: tmpstiff(ndof*nn, ndof*nn)
94 real(kind = kreal) :: elem(3, nn)
95 real(kind = kreal) :: unode(3, nn)
96 real(kind = kreal) :: xi_lx, eta_lx, zeta_ly
97 real(kind = kreal) :: w_w_lx, w_ly
98 real(kind = kreal) :: b_di(5, ndof*nn, 6, 3, 7)
99 real(kind = kreal) :: b1(3, ndof*nn), b2(3, ndof*nn), &
100 b3(3, ndof*nn)
101 real(kind = kreal) :: naturalcoord(2)
102 real(kind = kreal) :: tpcoord(6, 2, 3)
103 real(kind = kreal) :: nncoord(nn, 2)
104 real(kind = kreal) :: shapefunc(nn)
105 real(kind = kreal) :: shapederiv(nn, 2)
106 real(kind = kreal) :: aa1(3), aa2(3), aa3(3)
107 real(kind = kreal) :: bb1(3), bb2(3), bb3(3)
108 real(kind = kreal) :: cc1(3), cc2(3)
109 real(kind = kreal) :: alpha
110 real(kind = kreal) :: xxi_lx, eeta_lx
111 real(kind = kreal) :: xxi_di(6, 3), eeta_di(6, 3)
112 real(kind = kreal) :: h(nn, 3)
113 real(kind = kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
114 real(kind = kreal) :: v1_i(3), v2_i(3), v3_i(3)
115 real(kind = kreal) :: v1_abs, v2_abs, v3_abs
116 real(kind = kreal) :: a_over_2_v3(3, nn)
117 real(kind = kreal) :: u_rot(3, nn)
118 real(kind = kreal) :: dudxi_rot(3, nn), dudeta_rot(3, nn), &
119 dudzeta_rot(3, nn)
120 real(kind = kreal) :: g1(3), g2(3), g3(3)
121 real(kind = kreal) :: g3_abs
122 real(kind = kreal) :: e_0(3)
123 real(kind = kreal) :: cg1(3), cg2(3), cg3(3)
124 real(kind = kreal) :: det
125 real(kind = kreal) :: det_cg3(3)
126 real(kind = kreal) :: det_inv
127 real(kind = kreal) :: det_cg3_abs
128 real(kind = kreal) :: w_w_w_det
129 real(kind = kreal) :: e1_hat(3), e2_hat(3), e3_hat(3)
130 real(kind = kreal) :: e1_hat_abs, e2_hat_abs
131 real(kind = kreal) :: cv12(ndof*nn), cv13(ndof*nn), &
132 cv21(ndof*nn), cv23(ndof*nn), &
133 cv31(ndof*nn), cv32(ndof*nn)
134 real(kind = kreal) :: cv_theta(ndof*nn), cv_w(ndof*nn)
135 real(kind = kreal) :: cv(ndof*nn)
136 real(kind = kreal) :: sumlyr
137
138 sstable = 0
139 flag_dof = 0
140 ny = 0
141
142 !--------------------------------------------------------------------
143
144 ! MITC4
145 if( etype .EQ. fe_mitc4_shell ) then
146
147 fetype = fe_mitc4_shell
148
149 ny = 2
150
151 ntying = 1
152 npoints_tying(1)= 4
153
154 ! MITC9
155 else if( etype .EQ. fe_mitc9_shell ) then
156
157 fetype = fe_mitc9_shell
158
159 ny = 3
160
161 ntying = 3
162 npoints_tying(1)= 6
163 npoints_tying(2)= 6
164 npoints_tying(3)= 4
165
166 ! MITC3
167 else if( etype .EQ. fe_mitc3_shell ) then
168
169 fetype = fe_mitc3_shell
170
171 ny = 2
172
173 ntying = 1
174 npoints_tying(1)= 3
175
176 end if
177
178 !--------------------------------------------------------------------
179
180 if( present( nddisp ) ) then
181
182 unode(:, 1:nn) = nddisp(:, :)
183
184 end if
185
186 !--------------------------------------------------------------------
187
188 flag = gausses(1)%pMaterial%nlgeom_flag
189
190 if( .not. present( nddisp ) ) flag = infinitesimal
191
192 !--------------------------------------------------------------------
193
194 elem(:, :) = ecoord(:, :)
195
196 !--------------------------------------------------------------------
197
198 tmpstiff(:, :) = 0.0d0
199
200 !-------------------------------------------------------------------
201
202 ! xi-coordinate at a node in a local element
203 ! eta-coordinate at a node in a local element
204 call getnodalnaturalcoord(fetype, nncoord)
205
206 !-------------------------------------------------------------------
207
208 ! MITC4
209 if( etype .EQ. fe_mitc4_shell ) then
210
211 !--------------------------------------------------------
212
213 ! xi-coordinate at a tying pont in a local element
214 tpcoord(1, 1, 1) = 0.0d0
215 tpcoord(2, 1, 1) = 1.0d0
216 tpcoord(3, 1, 1) = 0.0d0
217 tpcoord(4, 1, 1) = -1.0d0
218 ! eta-coordinate at a tying point in a local element
219 tpcoord(1, 2, 1) = -1.0d0
220 tpcoord(2, 2, 1) = 0.0d0
221 tpcoord(3, 2, 1) = 1.0d0
222 tpcoord(4, 2, 1) = 0.0d0
223
224 !--------------------------------------------------------
225
226 ! MITC9
227 else if( etype .EQ. fe_mitc9_shell ) then
228
229 !--------------------------------------------------------
230
231 ! xi-coordinate at a tying point in a local element
232 tpcoord(1, 1, 1) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
233 tpcoord(2, 1, 1) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
234 tpcoord(3, 1, 1) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
235 tpcoord(4, 1, 1) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
236 tpcoord(5, 1, 1) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
237 tpcoord(6, 1, 1) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
238 ! eta-coordinate at a tying point in a local element
239 tpcoord(1, 2, 1) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
240 tpcoord(2, 2, 1) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
241 tpcoord(3, 2, 1) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
242 tpcoord(4, 2, 1) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
243 tpcoord(5, 2, 1) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
244 tpcoord(6, 2, 1) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
245
246 ! xi-coordinate at a tying point in a local element
247 tpcoord(1, 1, 2) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
248 tpcoord(2, 1, 2) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
249 tpcoord(3, 1, 2) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
250 tpcoord(4, 1, 2) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
251 tpcoord(5, 1, 2) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
252 tpcoord(6, 1, 2) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
253 ! eta-coordinate at a tying point in a local element
254 tpcoord(1, 2, 2) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
255 tpcoord(2, 2, 2) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
256 tpcoord(3, 2, 2) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
257 tpcoord(4, 2, 2) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
258 tpcoord(5, 2, 2) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
259 tpcoord(6, 2, 2) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
260
261 ! xi-coordinate at a tying point in a local element
262 tpcoord(1, 1, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
263 tpcoord(2, 1, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
264 tpcoord(3, 1, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
265 tpcoord(4, 1, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
266 ! eta-coordinate at a tying point in a local element
267 tpcoord(1, 2, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
268 tpcoord(2, 2, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
269 tpcoord(3, 2, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
270 tpcoord(4, 2, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
271
272 !--------------------------------------------------------
273
274 ! Xi-coordinate at a tying point in a local element
275 xxi_di(1, 1) = -1.0d0
276 xxi_di(2, 1) = 1.0d0
277 xxi_di(3, 1) = 1.0d0
278 xxi_di(4, 1) = -1.0d0
279 xxi_di(5, 1) = 1.0d0
280 xxi_di(6, 1) = -1.0d0
281 ! Eta-coordinate at a tying point in a local element
282 eeta_di(1, 1) = -1.0d0
283 eeta_di(2, 1) = -1.0d0
284 eeta_di(3, 1) = 1.0d0
285 eeta_di(4, 1) = 1.0d0
286 eeta_di(5, 1) = 0.0d0
287 eeta_di(6, 1) = 0.0d0
288
289 ! Xi-coordinate at a tying point in a local element
290 xxi_di(1, 2) = -1.0d0
291 xxi_di(2, 2) = 0.0d0
292 xxi_di(3, 2) = 1.0d0
293 xxi_di(4, 2) = 1.0d0
294 xxi_di(5, 2) = 0.0d0
295 xxi_di(6, 2) = -1.0d0
296 ! Eta-coordinate at a tying point in a local element
297 eeta_di(1, 2) = -1.0d0
298 eeta_di(2, 2) = -1.0d0
299 eeta_di(3, 2) = -1.0d0
300 eeta_di(4, 2) = 1.0d0
301 eeta_di(5, 2) = 1.0d0
302 eeta_di(6, 2) = 1.0d0
303
304 !--------------------------------------------------------
305
306 ! MITC3
307 else if( etype .EQ. fe_mitc3_shell ) then
308
309 !--------------------------------------------------------
310
311 ! xi-coordinate at a tying point in a local element
312 tpcoord(1, 1, 1) = 0.5d0
313 tpcoord(2, 1, 1) = 0.0d0
314 tpcoord(3, 1, 1) = 0.5d0
315 ! eta-coordinate at a tying point in a local element
316 tpcoord(1, 2, 1) = 0.0d0
317 tpcoord(2, 2, 1) = 0.5d0
318 tpcoord(3, 2, 1) = 0.5d0
319
320 !--------------------------------------------------------
321
322 end if
323
324 !--------------------------------------------------------------------
325
326 ! xi-coordinate at the center point in a local element
327 ! eta-coordinate at the center point in a local element
328 naturalcoord(1) = 0.0d0
329 naturalcoord(2) = 0.0d0
330
331 call getshapederiv(fetype, naturalcoord, shapederiv)
332
333 !--------------------------------------------------------------
334
335 ! Covariant basis vector
336 do i = 1, 3
337
338 g1(i) = 0.0d0
339
340 do na = 1, nn
341
342 g1(i) = g1(i)+shapederiv(na, 1) &
343 *elem(i, na)
344
345 end do
346
347 end do
348
349 e_0(1) = g1(1)
350 e_0(2) = g1(2)
351 e_0(3) = g1(3)
352
353 !--------------------------------------------------------------
354
355 do nb = 1, nn
356
357 !--------------------------------------------------------
358
359 naturalcoord(1) = nncoord(nb, 1)
360 naturalcoord(2) = nncoord(nb, 2)
361
362 call getshapederiv(fetype, naturalcoord, shapederiv)
363
364 !--------------------------------------------------------
365
366 ! Covariant basis vector
367 do i = 1, 3
368
369 g1(i) = 0.0d0
370 g2(i) = 0.0d0
371
372 do na = 1, nn
373
374 g1(i) = g1(i)+shapederiv(na, 1) &
375 *elem(i, na)
376 g2(i) = g2(i)+shapederiv(na, 2) &
377 *elem(i, na)
378
379 end do
380
381 end do
382
383 !------------------------------------------
384
385 det_cg3(1) = g1(2)*g2(3)-g1(3)*g2(2)
386 det_cg3(2) = g1(3)*g2(1)-g1(1)*g2(3)
387 det_cg3(3) = g1(1)*g2(2)-g1(2)*g2(1)
388
389 det_cg3_abs = dsqrt( det_cg3(1)*det_cg3(1) &
390 +det_cg3(2)*det_cg3(2) &
391 +det_cg3(3)*det_cg3(3) )
392
393 v3(1, nb) = det_cg3(1)/det_cg3_abs
394 v3(2, nb) = det_cg3(2)/det_cg3_abs
395 v3(3, nb) = det_cg3(3)/det_cg3_abs
396
397 !--------------------------------------------------------
398
399 v2(1, nb) = v3(2, nb)*e_0(3)-v3(3, nb)*e_0(2)
400 v2(2, nb) = v3(3, nb)*e_0(1)-v3(1, nb)*e_0(3)
401 v2(3, nb) = v3(1, nb)*e_0(2)-v3(2, nb)*e_0(1)
402
403 v2_abs = dsqrt( v2(1, nb)*v2(1, nb) &
404 +v2(2, nb)*v2(2, nb) &
405 +v2(3, nb)*v2(3, nb) )
406
407 if( v2_abs .GT. 1.0d-15 ) then
408
409 v2(1, nb) = v2(1, nb)/v2_abs
410 v2(2, nb) = v2(2, nb)/v2_abs
411 v2(3, nb) = v2(3, nb)/v2_abs
412
413 v1(1, nb) = v2(2, nb)*v3(3, nb) &
414 -v2(3, nb)*v3(2, nb)
415 v1(2, nb) = v2(3, nb)*v3(1, nb) &
416 -v2(1, nb)*v3(3, nb)
417 v1(3, nb) = v2(1, nb)*v3(2, nb) &
418 -v2(2, nb)*v3(1, nb)
419
420 v1_abs = dsqrt( v1(1, nb)*v1(1, nb) &
421 +v1(2, nb)*v1(2, nb) &
422 +v1(3, nb)*v1(3, nb) )
423
424 v1(1, nb) = v1(1, nb)/v1_abs
425 v1(2, nb) = v1(2, nb)/v1_abs
426 v1(3, nb) = v1(3, nb)/v1_abs
427
428 else ! YX: impossible
429
430 v1(1, nb) = 0.0d0
431 v1(2, nb) = 0.0d0
432 v1(3, nb) = -1.0d0
433
434 v2(1, nb) = 0.0d0
435 v2(2, nb) = 1.0d0
436 v2(3, nb) = 0.0d0
437
438 end if
439
440 !---------------------------------------------------
441
442 v3(1, nb) = v1(2, nb)*v2(3, nb) &
443 -v1(3, nb)*v2(2, nb)
444 v3(2, nb) = v1(3, nb)*v2(1, nb) &
445 -v1(1, nb)*v2(3, nb)
446 v3(3, nb) = v1(1, nb)*v2(2, nb) &
447 -v1(2, nb)*v2(1, nb)
448
449 v3_abs = dsqrt( v3(1, nb)*v3(1, nb) &
450 +v3(2, nb)*v3(2, nb) &
451 +v3(3, nb)*v3(3, nb) )
452
453 v3(1, nb) = v3(1, nb)/v3_abs
454 v3(2, nb) = v3(2, nb)/v3_abs
455 v3(3, nb) = v3(3, nb)/v3_abs
456
457 !--------------------------------------------------------
458
459 a_over_2_v3(1, nb) = 0.5d0*thick*v3(1, nb)
460 a_over_2_v3(2, nb) = 0.5d0*thick*v3(2, nb)
461 a_over_2_v3(3, nb) = 0.5d0*thick*v3(3, nb)
462
463 !--------------------------------------------------------
464
465 end do
466
467 !--------------------------------------------------------------------
468 ! MODIFIED to LAMINATED SHELL ANALYSIS
469 !--------------------------------------------------------------------
470 zeta_ly = 0.0d0
471
472 n_totlyr = gausses(1)%pMaterial%totallyr
473 do n_layer=1,n_totlyr
474 do ly = 1, ny
475
476 !--------------------------------------------------------
477
478 ! MITC4
479 if( etype .EQ. fe_mitc4_shell ) then
480
481 zeta_ly = 0.0d0
482
483 ! MITC9
484 else if( etype .EQ. fe_mitc9_shell ) then
485
486 zeta_ly = xg(ny, ly)
487
488 ! MITC3
489 else if( etype .EQ. fe_mitc3_shell )then
490
491 zeta_ly = 0.0d0
492
493 end if
494
495 !--------------------------------------------------------
496
497 do it = 1, ntying
498
499 do ip = 1, npoints_tying(it)
500
501 !-------------------------------------------------
502
503 naturalcoord(1) = tpcoord(ip, 1, it)
504 naturalcoord(2) = tpcoord(ip, 2, it)
505
506 call getshapefunc(fetype, naturalcoord, shapefunc)
507
508 call getshapederiv(fetype, naturalcoord, shapederiv)
509
510 !-------------------------------------------------
511
512 do na = 1, nn
513
514 do i = 1, 3
515
516 u_rot(i, na) &
517 = shapefunc(na) &
518 *( zeta_ly*a_over_2_v3(i, na) )
519
520 dudxi_rot(i, na) &
521 = shapederiv(na, 1) &
522 *( zeta_ly*a_over_2_v3(i, na) )
523 dudeta_rot(i, na) &
524 = shapederiv(na, 2) &
525 *( zeta_ly*a_over_2_v3(i, na) )
526 dudzeta_rot(i, na) &
527 = shapefunc(na) &
528 *( a_over_2_v3(i, na) )
529
530 end do
531
532 end do
533
534 !-------------------------------------------------
535
536 ! Covariant basis vector
537 do i = 1, 3
538
539 g1(i) = 0.0d0
540 g2(i) = 0.0d0
541 g3(i) = 0.0d0
542
543 do na = 1, nn
544
545 g1(i) = g1(i)+shapederiv(na, 1) &
546 *elem(i, na) &
547 +dudxi_rot(i, na)
548 g2(i) = g2(i)+shapederiv(na, 2) &
549 *elem(i, na) &
550 +dudeta_rot(i, na)
551 g3(i) = g3(i)+dudzeta_rot(i, na)
552
553 end do
554
555 end do
556
557 !-------------------------------------------------
558
559 ! [ B L ] matrix
560 do nb = 1, nn
561
562 jsize1 = ndof*(nb-1)+1
563 jsize2 = ndof*(nb-1)+2
564 jsize3 = ndof*(nb-1)+3
565 jsize4 = ndof*(nb-1)+4
566 jsize5 = ndof*(nb-1)+5
567 jsize6 = ndof*(nb-1)+6
568
569 aa1(1) = dudxi_rot(2, nb) *g1(3)-dudxi_rot(3, nb) *g1(2)
570 aa1(2) = dudxi_rot(3, nb) *g1(1)-dudxi_rot(1, nb) *g1(3)
571 aa1(3) = dudxi_rot(1, nb) *g1(2)-dudxi_rot(2, nb) *g1(1)
572
573 aa2(1) = dudxi_rot(2, nb) *g2(3)-dudxi_rot(3, nb) *g2(2)
574 aa2(2) = dudxi_rot(3, nb) *g2(1)-dudxi_rot(1, nb) *g2(3)
575 aa2(3) = dudxi_rot(1, nb) *g2(2)-dudxi_rot(2, nb) *g2(1)
576
577 aa3(1) = dudxi_rot(2, nb) *g3(3)-dudxi_rot(3, nb) *g3(2)
578 aa3(2) = dudxi_rot(3, nb) *g3(1)-dudxi_rot(1, nb) *g3(3)
579 aa3(3) = dudxi_rot(1, nb) *g3(2)-dudxi_rot(2, nb) *g3(1)
580
581 bb1(1) = dudeta_rot(2, nb) *g1(3)-dudeta_rot(3, nb) *g1(2)
582 bb1(2) = dudeta_rot(3, nb) *g1(1)-dudeta_rot(1, nb) *g1(3)
583 bb1(3) = dudeta_rot(1, nb) *g1(2)-dudeta_rot(2, nb) *g1(1)
584
585 bb2(1) = dudeta_rot(2, nb) *g2(3)-dudeta_rot(3, nb) *g2(2)
586 bb2(2) = dudeta_rot(3, nb) *g2(1)-dudeta_rot(1, nb) *g2(3)
587 bb2(3) = dudeta_rot(1, nb) *g2(2)-dudeta_rot(2, nb) *g2(1)
588
589 bb3(1) = dudeta_rot(2, nb) *g3(3)-dudeta_rot(3, nb) *g3(2)
590 bb3(2) = dudeta_rot(3, nb) *g3(1)-dudeta_rot(1, nb) *g3(3)
591 bb3(3) = dudeta_rot(1, nb) *g3(2)-dudeta_rot(2, nb) *g3(1)
592
593 cc1(1) = dudzeta_rot(2, nb)*g1(3)-dudzeta_rot(3, nb)*g1(2)
594 cc1(2) = dudzeta_rot(3, nb)*g1(1)-dudzeta_rot(1, nb)*g1(3)
595 cc1(3) = dudzeta_rot(1, nb)*g1(2)-dudzeta_rot(2, nb)*g1(1)
596
597 cc2(1) = dudzeta_rot(2, nb)*g2(3)-dudzeta_rot(3, nb)*g2(2)
598 cc2(2) = dudzeta_rot(3, nb)*g2(1)-dudzeta_rot(1, nb)*g2(3)
599 cc2(3) = dudzeta_rot(1, nb)*g2(2)-dudzeta_rot(2, nb)*g2(1)
600
601 b_di(1, jsize1, ip, it, ly) = shapederiv(nb, 1)*g1(1)
602 b_di(2, jsize1, ip, it, ly) = shapederiv(nb, 2)*g2(1)
603 b_di(3, jsize1, ip, it, ly) = shapederiv(nb, 1)*g2(1) &
604 +shapederiv(nb, 2)*g1(1)
605 b_di(4, jsize1, ip, it, ly) = shapederiv(nb, 2)*g3(1)
606 b_di(5, jsize1, ip, it, ly) = shapederiv(nb, 1)*g3(1)
607
608 b_di(1, jsize2, ip, it, ly) = shapederiv(nb, 1)*g1(2)
609 b_di(2, jsize2, ip, it, ly) = shapederiv(nb, 2)*g2(2)
610 b_di(3, jsize2, ip, it, ly) = shapederiv(nb, 1)*g2(2) &
611 +shapederiv(nb, 2)*g1(2)
612 b_di(4, jsize2, ip, it, ly) = shapederiv(nb, 2)*g3(2)
613 b_di(5, jsize2, ip, it, ly) = shapederiv(nb, 1)*g3(2)
614
615 b_di(1, jsize3, ip, it, ly) = shapederiv(nb, 1)*g1(3)
616 b_di(2, jsize3, ip, it, ly) = shapederiv(nb, 2)*g2(3)
617 b_di(3, jsize3, ip, it, ly) = shapederiv(nb, 1)*g2(3) &
618 +shapederiv(nb, 2)*g1(3)
619 b_di(4, jsize3, ip, it, ly) = shapederiv(nb, 2)*g3(3)
620 b_di(5, jsize3, ip, it, ly) = shapederiv(nb, 1)*g3(3)
621
622 b_di(1, jsize4, ip, it, ly) = aa1(1)
623 b_di(2, jsize4, ip, it, ly) = bb2(1)
624 b_di(3, jsize4, ip, it, ly) = aa2(1)+bb1(1)
625 b_di(4, jsize4, ip, it, ly) = bb3(1)+cc2(1)
626 b_di(5, jsize4, ip, it, ly) = aa3(1)+cc1(1)
627
628 b_di(1, jsize5, ip, it, ly) = aa1(2)
629 b_di(2, jsize5, ip, it, ly) = bb2(2)
630 b_di(3, jsize5, ip, it, ly) = aa2(2)+bb1(2)
631 b_di(4, jsize5, ip, it, ly) = bb3(2)+cc2(2)
632 b_di(5, jsize5, ip, it, ly) = aa3(2)+cc1(2)
633
634 b_di(1, jsize6, ip, it, ly) = aa1(3)
635 b_di(2, jsize6, ip, it, ly) = bb2(3)
636 b_di(3, jsize6, ip, it, ly) = aa2(3)+bb1(3)
637 b_di(4, jsize6, ip, it, ly) = bb3(3)+cc2(3)
638 b_di(5, jsize6, ip, it, ly) = aa3(3)+cc1(3)
639
640 end do
641
642 !-------------------------------------------------
643
644 end do
645
646 end do
647
648 !--------------------------------------------------------
649
650 sumlyr = 0.0d0
651 do m = 1, n_layer
652 sumlyr = sumlyr +2 *gausses(1)%pMaterial%shell_var(m)%weight
653 enddo
654 zeta_ly = -1 +sumlyr -gausses(1)%pMaterial%shell_var(n_layer)%weight*(1-xg(ny, ly))
655
656 w_ly = wgt(ny, ly)
657
658 !--------------------------------------------------------
659
660 do lx = 1, numofquadpoints(fetype)
661
662 !--------------------------------------------------
663
664 call getquadpoint(fetype, lx, naturalcoord)
665
666 xi_lx = naturalcoord(1)
667 eta_lx = naturalcoord(2)
668
669 w_w_lx = getweight(fetype, lx)
670
671 call getshapefunc(fetype, naturalcoord, shapefunc)
672
673 call getshapederiv(fetype, naturalcoord, shapederiv)
674
675 !--------------------------------------------------
676
677 do i = 1, 3
678
679 v1_i(i) = 0.0d0
680 v2_i(i) = 0.0d0
681 v3_i(i) = 0.0d0
682
683 do na = 1, nn
684
685 v1_i(i) = v1_i(i)+shapefunc(na)*v1(i, na)
686 v2_i(i) = v2_i(i)+shapefunc(na)*v2(i, na)
687 v3_i(i) = v3_i(i)+shapefunc(na)*v3(i, na)
688
689 end do
690
691 end do
692
693 !--------------------------------------------------
694
695 do na = 1, nn
696
697 do i = 1, 3
698
699 u_rot(i, na) &
700 = shapefunc(na) &
701 *( zeta_ly*a_over_2_v3(i, na) )
702
703 dudxi_rot(i, na) &
704 = shapederiv(na, 1) &
705 *( zeta_ly*a_over_2_v3(i, na) )
706 dudeta_rot(i, na) &
707 = shapederiv(na, 2) &
708 *( zeta_ly*a_over_2_v3(i, na) )
709 dudzeta_rot(i, na) &
710 = shapefunc(na) &
711 *( a_over_2_v3(i, na) )
712
713 end do
714
715 end do
716
717 !--------------------------------------------------
718
719 ! Covariant basis vector
720 do i = 1, 3
721
722 g1(i) = 0.0d0
723 g2(i) = 0.0d0
724 g3(i) = 0.0d0
725
726 do na = 1, nn
727
728 g1(i) = g1(i)+shapederiv(na, 1) &
729 *elem(i, na) &
730 +dudxi_rot(i, na)
731 g2(i) = g2(i)+shapederiv(na, 2) &
732 *elem(i, na) &
733 +dudeta_rot(i, na)
734 g3(i) = g3(i)+dudzeta_rot(i, na)
735
736 end do
737
738 end do
739
740 !--------------------------------------------------
741
742 ! Jacobian
743 det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
744 +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
745 +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )
746
747 det_inv = 1.0d0/det
748
749 !--------------------------------------------------
750
751 ! Contravariant basis vector
752 cg1(1) = det_inv &
753 *( g2(2)*g3(3)-g2(3)*g3(2) )
754 cg1(2) = det_inv &
755 *( g2(3)*g3(1)-g2(1)*g3(3) )
756 cg1(3) = det_inv &
757 *( g2(1)*g3(2)-g2(2)*g3(1) )
758 cg2(1) = det_inv &
759 *( g3(2)*g1(3)-g3(3)*g1(2) )
760 cg2(2) = det_inv &
761 *( g3(3)*g1(1)-g3(1)*g1(3) )
762 cg2(3) = det_inv &
763 *( g3(1)*g1(2)-g3(2)*g1(1) )
764 cg3(1) = det_inv &
765 *( g1(2)*g2(3)-g1(3)*g2(2) )
766 cg3(2) = det_inv &
767 *( g1(3)*g2(1)-g1(1)*g2(3) )
768 cg3(3) = det_inv &
769 *( g1(1)*g2(2)-g1(2)*g2(1) )
770
771 !--------------------------------------------------
772
773 g3_abs = dsqrt( g3(1)*g3(1) &
774 +g3(2)*g3(2) &
775 +g3(3)*g3(3) )
776
777 !--------------------------------------------------
778
779 ! Orthonormal vectors
780
781 e3_hat(1) = g3(1)/g3_abs
782 e3_hat(2) = g3(2)/g3_abs
783 e3_hat(3) = g3(3)/g3_abs
784
785 e1_hat(1) = g2(2)*e3_hat(3) &
786 -g2(3)*e3_hat(2)
787 e1_hat(2) = g2(3)*e3_hat(1) &
788 -g2(1)*e3_hat(3)
789 e1_hat(3) = g2(1)*e3_hat(2) &
790 -g2(2)*e3_hat(1)
791 e1_hat_abs = dsqrt( e1_hat(1)*e1_hat(1) &
792 +e1_hat(2)*e1_hat(2) &
793 +e1_hat(3)*e1_hat(3) )
794 e1_hat(1) = e1_hat(1)/e1_hat_abs
795 e1_hat(2) = e1_hat(2)/e1_hat_abs
796 e1_hat(3) = e1_hat(3)/e1_hat_abs
797
798 e2_hat(1) = e3_hat(2)*e1_hat(3) &
799 -e3_hat(3)*e1_hat(2)
800 e2_hat(2) = e3_hat(3)*e1_hat(1) &
801 -e3_hat(1)*e1_hat(3)
802 e2_hat(3) = e3_hat(1)*e1_hat(2) &
803 -e3_hat(2)*e1_hat(1)
804 e2_hat_abs = dsqrt( e2_hat(1)*e2_hat(1) &
805 +e2_hat(2)*e2_hat(2) &
806 +e2_hat(3)*e2_hat(3) )
807 e2_hat(1) = e2_hat(1)/e2_hat_abs
808 e2_hat(2) = e2_hat(2)/e2_hat_abs
809 e2_hat(3) = e2_hat(3)/e2_hat_abs
810
811 !--------------------------------------------------
812
813 call matlmatrix_shell &
814 (gausses(lx), shell, d, &
815 e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
816 alpha, n_layer)
817
818 !--------------------------------------------------
819
820 ! [ B L ] matrix
821 do nb = 1, nn
822
823 jsize1 = ndof*(nb-1)+1
824 jsize2 = ndof*(nb-1)+2
825 jsize3 = ndof*(nb-1)+3
826 jsize4 = ndof*(nb-1)+4
827 jsize5 = ndof*(nb-1)+5
828 jsize6 = ndof*(nb-1)+6
829
830 aa1(1) = dudxi_rot(2, nb) *g1(3)-dudxi_rot(3, nb) *g1(2)
831 aa1(2) = dudxi_rot(3, nb) *g1(1)-dudxi_rot(1, nb) *g1(3)
832 aa1(3) = dudxi_rot(1, nb) *g1(2)-dudxi_rot(2, nb) *g1(1)
833
834 aa2(1) = dudxi_rot(2, nb) *g2(3)-dudxi_rot(3, nb) *g2(2)
835 aa2(2) = dudxi_rot(3, nb) *g2(1)-dudxi_rot(1, nb) *g2(3)
836 aa2(3) = dudxi_rot(1, nb) *g2(2)-dudxi_rot(2, nb) *g2(1)
837
838 aa3(1) = dudxi_rot(2, nb) *g3(3)-dudxi_rot(3, nb) *g3(2)
839 aa3(2) = dudxi_rot(3, nb) *g3(1)-dudxi_rot(1, nb) *g3(3)
840 aa3(3) = dudxi_rot(1, nb) *g3(2)-dudxi_rot(2, nb) *g3(1)
841
842 bb1(1) = dudeta_rot(2, nb) *g1(3)-dudeta_rot(3, nb) *g1(2)
843 bb1(2) = dudeta_rot(3, nb) *g1(1)-dudeta_rot(1, nb) *g1(3)
844 bb1(3) = dudeta_rot(1, nb) *g1(2)-dudeta_rot(2, nb) *g1(1)
845
846 bb2(1) = dudeta_rot(2, nb) *g2(3)-dudeta_rot(3, nb) *g2(2)
847 bb2(2) = dudeta_rot(3, nb) *g2(1)-dudeta_rot(1, nb) *g2(3)
848 bb2(3) = dudeta_rot(1, nb) *g2(2)-dudeta_rot(2, nb) *g2(1)
849
850 bb3(1) = dudeta_rot(2, nb) *g3(3)-dudeta_rot(3, nb) *g3(2)
851 bb3(2) = dudeta_rot(3, nb) *g3(1)-dudeta_rot(1, nb) *g3(3)
852 bb3(3) = dudeta_rot(1, nb) *g3(2)-dudeta_rot(2, nb) *g3(1)
853
854 cc1(1) = dudzeta_rot(2, nb)*g1(3)-dudzeta_rot(3, nb)*g1(2)
855 cc1(2) = dudzeta_rot(3, nb)*g1(1)-dudzeta_rot(1, nb)*g1(3)
856 cc1(3) = dudzeta_rot(1, nb)*g1(2)-dudzeta_rot(2, nb)*g1(1)
857
858 cc2(1) = dudzeta_rot(2, nb)*g2(3)-dudzeta_rot(3, nb)*g2(2)
859 cc2(2) = dudzeta_rot(3, nb)*g2(1)-dudzeta_rot(1, nb)*g2(3)
860 cc2(3) = dudzeta_rot(1, nb)*g2(2)-dudzeta_rot(2, nb)*g2(1)
861
862 b(1, jsize1) = shapederiv(nb, 1)*g1(1)
863 b(2, jsize1) = shapederiv(nb, 2)*g2(1)
864 b(3, jsize1) = shapederiv(nb, 1)*g2(1) &
865 +shapederiv(nb, 2)*g1(1)
866 b(4, jsize1) = shapederiv(nb, 2)*g3(1)
867 b(5, jsize1) = shapederiv(nb, 1)*g3(1)
868
869 b(1, jsize2) = shapederiv(nb, 1)*g1(2)
870 b(2, jsize2) = shapederiv(nb, 2)*g2(2)
871 b(3, jsize2) = shapederiv(nb, 1)*g2(2) &
872 +shapederiv(nb, 2)*g1(2)
873 b(4, jsize2) = shapederiv(nb, 2)*g3(2)
874 b(5, jsize2) = shapederiv(nb, 1)*g3(2)
875
876 b(1, jsize3) = shapederiv(nb, 1)*g1(3)
877 b(2, jsize3) = shapederiv(nb, 2)*g2(3)
878 b(3, jsize3) = shapederiv(nb, 1)*g2(3) &
879 +shapederiv(nb, 2)*g1(3)
880 b(4, jsize3) = shapederiv(nb, 2)*g3(3)
881 b(5, jsize3) = shapederiv(nb, 1)*g3(3)
882
883 b(1, jsize4) = aa1(1)
884 b(2, jsize4) = bb2(1)
885 b(3, jsize4) = aa2(1)+bb1(1)
886 b(4, jsize4) = bb3(1)+cc2(1)
887 b(5, jsize4) = aa3(1)+cc1(1)
888
889 b(1, jsize5) = aa1(2)
890 b(2, jsize5) = bb2(2)
891 b(3, jsize5) = aa2(2)+bb1(2)
892 b(4, jsize5) = bb3(2)+cc2(2)
893 b(5, jsize5) = aa3(2)+cc1(2)
894
895 b(1, jsize6) = aa1(3)
896 b(2, jsize6) = bb2(3)
897 b(3, jsize6) = aa2(3)+bb1(3)
898 b(4, jsize6) = bb3(3)+cc2(3)
899 b(5, jsize6) = aa3(3)+cc1(3)
900
901 end do
902
903 !--------------------------------------------------
904
905 ! MITC4
906 if( etype .EQ. fe_mitc4_shell ) then
907
908 do jsize = 1, ndof*nn
909
910 b(4, jsize) = 0.0d0
911 b(5, jsize) = 0.0d0
912
913 ! B_as(4, jsize)
914 b(4, jsize) &
915 = 0.5d0*( 1.0d0-xi_lx )*b_di(4, jsize, 4, 1, ly) &
916 +0.5d0*( 1.0d0+xi_lx )*b_di(4, jsize, 2, 1, ly)
917 ! B_as(5, jsize)
918 b(5, jsize) &
919 = 0.5d0*( 1.0d0-eta_lx )*b_di(5, jsize, 1, 1, ly) &
920 +0.5d0*( 1.0d0+eta_lx )*b_di(5, jsize, 3, 1, ly)
921
922 end do
923
924 ! MITC9
925 else if( etype .EQ. fe_mitc9_shell ) then
926
927 xxi_lx = xi_lx /dsqrt( 1.0d0/3.0d0 )
928 eeta_lx = eta_lx/dsqrt( 3.0d0/5.0d0 )
929
930 do ip = 1, npoints_tying(1)
931
932 h(ip, 1) &
933 = ( 0.5d0*( 1.0d0+xxi_di(ip, 1)*xxi_lx ) ) &
934 *( ( 0.5d0*eeta_di(ip, 1)*eeta_lx ) &
935 *( 1.0d0+eeta_di(ip, 1)*eeta_lx ) &
936 +( 1.0d0-eeta_di(ip, 1)*eeta_di(ip, 1) ) &
937 *( 1.0d0-eeta_lx*eeta_lx ) )
938
939 end do
940
941 xxi_lx = xi_lx /dsqrt( 3.0d0/5.0d0 )
942 eeta_lx = eta_lx/dsqrt( 1.0d0/3.0d0 )
943
944 do ip = 1, npoints_tying(2)
945
946 h(ip, 2) &
947 = ( ( 0.5d0*xxi_di(ip, 2) *xxi_lx ) &
948 *( 1.0d0+xxi_di(ip, 2) *xxi_lx ) &
949 +( 1.0d0-xxi_di(ip, 2) *xxi_di(ip, 2) ) &
950 *( 1.0d0-xxi_lx*xxi_lx ) ) &
951 *( 0.5d0*( 1.0d0+eeta_di(ip, 2)*eeta_lx ) )
952
953 end do
954
955 xxi_lx = xi_lx /dsqrt( 1.0d0/3.0d0 )
956 eeta_lx = eta_lx/dsqrt( 1.0d0/3.0d0 )
957
958 do ip = 1, npoints_tying(3)
959
960 h(ip, 3) &
961 = ( 0.5d0*( 1.0d0+xxi_di(ip, 1)*xxi_lx ) ) &
962 *( 0.5d0*( 1.0d0+eeta_di(ip, 1)*eeta_lx ) )
963
964 end do
965
966 do jsize = 1, ndof*nn
967
968 b(1, jsize) = 0.0d0
969 b(2, jsize) = 0.0d0
970 b(3, jsize) = 0.0d0
971 b(4, jsize) = 0.0d0
972 b(5, jsize) = 0.0d0
973
974 do ip = 1, npoints_tying(1)
975
976 ! B_as(1, jsize)
977 b(1, jsize) &
978 = b(1, jsize)+h(ip, 1)*b_di(1, jsize, ip, 1, ly)
979 ! B_as(5, jsize)
980 b(5, jsize) &
981 = b(5, jsize)+h(ip, 1)*b_di(5, jsize, ip, 1, ly)
982
983 end do
984
985 do ip = 1, npoints_tying(2)
986
987 ! B_as(2, jsize)
988 b(2, jsize) &
989 = b(2, jsize)+h(ip, 2)*b_di(2, jsize, ip, 2, ly)
990 ! B_as(4, jsize)
991 b(4, jsize) &
992 = b(4, jsize)+h(ip, 2)*b_di(4, jsize, ip, 2, ly)
993
994 end do
995
996 do ip = 1, npoints_tying(3)
997
998 ! B_as(3, jsize)
999 b(3, jsize) &
1000 = b(3, jsize)+h(ip, 3)*b_di(3, jsize, ip, 3, ly)
1001
1002 end do
1003
1004 end do
1005
1006 ! MITC3
1007 else if( etype .EQ. fe_mitc3_shell ) then
1008
1009 do jsize = 1, ndof*nn
1010
1011 b(4, jsize) = 0.0d0
1012 b(5, jsize) = 0.0d0
1013
1014 ! B_as(4, jsize)
1015 b(4, jsize) &
1016 = ( 1.0d0-xi_lx )*b_di(4, jsize, 2, 1, ly) &
1017 +xi_lx *b_di(5, jsize, 1, 1, ly) &
1018 +xi_lx *( b_di(4, jsize, 3, 1, ly) &
1019 -b_di(5, jsize, 3, 1, ly) )
1020
1021 ! B_as(5, jsize)
1022 b(5, jsize) &
1023 = eta_lx*b_di(4, jsize, 2, 1, ly) &
1024 +( 1.0d0-eta_lx )*b_di(5, jsize, 1, 1, ly) &
1025 -eta_lx*( b_di(4, jsize, 3, 1, ly) &
1026 -b_di(5, jsize, 3, 1, ly) )
1027
1028 end do
1029
1030 end if
1031
1032 !--------------------------------------------------
1033
1034 w_w_w_det = w_w_lx*w_ly*det
1035
1036 !--------------------------------------------------
1037
1038 db(1:5, 1:ndof*nn) = matmul( d, b(1:5, 1:ndof*nn ) )
1039
1040 !--------------------------------------------------
1041
1042 forall( isize=1:ndof*nn, jsize=1:ndof*nn )
1043 tmpstiff(isize, jsize) &
1044 = tmpstiff(isize, jsize) &
1045 +w_w_w_det*gausses(1)%pMaterial%shell_var(n_layer)%weight*dot_product( b(:, isize), db(:, jsize) )
1046 end forall
1047
1048 !--------------------------------------------------
1049
1050 ! [ B_{i} ] matrix
1051 do nb = 1, nn
1052
1053 jsize1 = ndof*(nb-1)+1
1054 jsize2 = ndof*(nb-1)+2
1055 jsize3 = ndof*(nb-1)+3
1056 jsize4 = ndof*(nb-1)+4
1057 jsize5 = ndof*(nb-1)+5
1058 jsize6 = ndof*(nb-1)+6
1059
1060 b1(1, jsize1) = shapederiv(nb, 1)
1061 b1(2, jsize1) = 0.0d0
1062 b1(3, jsize1) = 0.0d0
1063 b1(1, jsize2) = 0.0d0
1064 b1(2, jsize2) = shapederiv(nb, 1)
1065 b1(3, jsize2) = 0.0d0
1066 b1(1, jsize3) = 0.0d0
1067 b1(2, jsize3) = 0.0d0
1068 b1(3, jsize3) = shapederiv(nb, 1)
1069 b1(1, jsize4) = 0.0d0
1070 b1(2, jsize4) = -dudxi_rot(3, nb)
1071 b1(3, jsize4) = dudxi_rot(2, nb)
1072 b1(1, jsize5) = dudxi_rot(3, nb)
1073 b1(2, jsize5) = 0.0d0
1074 b1(3, jsize5) = -dudxi_rot(1, nb)
1075 b1(1, jsize6) = -dudxi_rot(2, nb)
1076 b1(2, jsize6) = dudxi_rot(1, nb)
1077 b1(3, jsize6) = 0.0d0
1078
1079 b2(1, jsize1) = shapederiv(nb, 2)
1080 b2(2, jsize1) = 0.0d0
1081 b2(3, jsize1) = 0.0d0
1082 b2(1, jsize2) = 0.0d0
1083 b2(2, jsize2) = shapederiv(nb, 2)
1084 b2(3, jsize2) = 0.0d0
1085 b2(1, jsize3) = 0.0d0
1086 b2(2, jsize3) = 0.0d0
1087 b2(3, jsize3) = shapederiv(nb, 2)
1088 b2(1, jsize4) = 0.0d0
1089 b2(2, jsize4) = -dudeta_rot(3, nb)
1090 b2(3, jsize4) = dudeta_rot(2, nb)
1091 b2(1, jsize5) = dudeta_rot(3, nb)
1092 b2(2, jsize5) = 0.0d0
1093 b2(3, jsize5) = -dudeta_rot(1, nb)
1094 b2(1, jsize6) = -dudeta_rot(2, nb)
1095 b2(2, jsize6) = dudeta_rot(1, nb)
1096 b2(3, jsize6) = 0.0d0
1097
1098 b3(1, jsize1) = 0.0d0
1099 b3(2, jsize1) = 0.0d0
1100 b3(3, jsize1) = 0.0d0
1101 b3(1, jsize2) = 0.0d0
1102 b3(2, jsize2) = 0.0d0
1103 b3(3, jsize2) = 0.0d0
1104 b3(1, jsize3) = 0.0d0
1105 b3(2, jsize3) = 0.0d0
1106 b3(3, jsize3) = 0.0d0
1107 b3(1, jsize4) = 0.0d0
1108 b3(2, jsize4) = -dudzeta_rot(3, nb)
1109 b3(3, jsize4) = dudzeta_rot(2, nb)
1110 b3(1, jsize5) = dudzeta_rot(3, nb)
1111 b3(2, jsize5) = 0.0d0
1112 b3(3, jsize5) = -dudzeta_rot(1, nb)
1113 b3(1, jsize6) = -dudzeta_rot(2, nb)
1114 b3(2, jsize6) = dudzeta_rot(1, nb)
1115 b3(3, jsize6) = 0.0d0
1116
1117 end do
1118
1119 !--------------------------------------------------
1120
1121 ! { C_{ij} } vector
1122 do jsize = 1, ndof*nn
1123
1124 cv12(jsize) = ( cg1(1)*b1(2, jsize) &
1125 +cg2(1)*b2(2, jsize) &
1126 +cg3(1)*b3(2, jsize) ) &
1127 -( cg1(2)*b1(1, jsize) &
1128 +cg2(2)*b2(1, jsize) &
1129 +cg3(2)*b3(1, jsize) )
1130 cv13(jsize) = ( cg1(1)*b1(3, jsize) &
1131 +cg2(1)*b2(3, jsize) &
1132 +cg3(1)*b3(3, jsize) ) &
1133 -( cg1(3)*b1(1, jsize) &
1134 +cg2(3)*b2(1, jsize) &
1135 +cg3(3)*b3(1, jsize) )
1136 cv21(jsize) = ( cg1(2)*b1(1, jsize) &
1137 +cg2(2)*b2(1, jsize) &
1138 +cg3(2)*b3(1, jsize) ) &
1139 -( cg1(1)*b1(2, jsize) &
1140 +cg2(1)*b2(2, jsize) &
1141 +cg3(1)*b3(2, jsize) )
1142 cv23(jsize) = ( cg1(2)*b1(3, jsize) &
1143 +cg2(2)*b2(3, jsize) &
1144 +cg3(2)*b3(3, jsize) ) &
1145 -( cg1(3)*b1(2, jsize) &
1146 +cg2(3)*b2(2, jsize) &
1147 +cg3(3)*b3(2, jsize) )
1148 cv31(jsize) = ( cg1(3)*b1(1, jsize) &
1149 +cg2(3)*b2(1, jsize) &
1150 +cg3(3)*b3(1, jsize) ) &
1151 -( cg1(1)*b1(3, jsize) &
1152 +cg2(1)*b2(3, jsize) &
1153 +cg3(1)*b3(3, jsize) )
1154 cv32(jsize) = ( cg1(3)*b1(2, jsize) &
1155 +cg2(3)*b2(2, jsize) &
1156 +cg3(3)*b3(2, jsize) ) &
1157 -( cg1(2)*b1(3, jsize) &
1158 +cg2(2)*b2(3, jsize) &
1159 +cg3(2)*b3(3, jsize) )
1160
1161 end do
1162
1163 !--------------------------------------------------
1164
1165 ! { Cw } vector
1166 do nb = 1, nn
1167
1168 do j = 1, ndof
1169
1170 jsize = ndof*(nb-1)+j
1171
1172 cv_w(jsize) &
1173 = v1_i(1)*cv12(jsize)*v2_i(2) &
1174 +v1_i(1)*cv13(jsize)*v2_i(3) &
1175 +v1_i(2)*cv21(jsize)*v2_i(1) &
1176 +v1_i(2)*cv23(jsize)*v2_i(3) &
1177 +v1_i(3)*cv31(jsize)*v2_i(1) &
1178 +v1_i(3)*cv32(jsize)*v2_i(2)
1179
1180 end do
1181
1182 end do
1183
1184 ! { Ctheta } vector
1185 do nb = 1, nn
1186
1187 jsize1 = ndof*(nb-1)+1
1188 jsize2 = ndof*(nb-1)+2
1189 jsize3 = ndof*(nb-1)+3
1190 jsize4 = ndof*(nb-1)+4
1191 jsize5 = ndof*(nb-1)+5
1192 jsize6 = ndof*(nb-1)+6
1193
1194 cv_theta(jsize1) = 0.0d0
1195 cv_theta(jsize2) = 0.0d0
1196 cv_theta(jsize3) = 0.0d0
1197 cv_theta(jsize4) = v3_i(1)*shapefunc(nb)
1198 cv_theta(jsize5) = v3_i(2)*shapefunc(nb)
1199 cv_theta(jsize6) = v3_i(3)*shapefunc(nb)
1200
1201 end do
1202
1203 ! { C } vector
1204 do jsize = 1, ndof*nn
1205
1206 cv(jsize) = cv_theta(jsize)-0.5d0*cv_w(jsize)
1207
1208 end do
1209
1210 !--------------------------------------------------
1211
1212 ! [ K L ] matrix
1213 do jsize = 1, ndof*nn
1214 do isize = 1, ndof*nn
1215
1216 tmpstiff(isize, jsize) &
1217 = tmpstiff(isize, jsize) &
1218 +w_w_w_det*gausses(1)%pMaterial%shell_var(n_layer)%weight*alpha*cv(isize)*cv(jsize)
1219
1220 end do
1221 end do
1222
1223
1224 !--------------------------------------------------
1225
1226 end do
1227
1228 !--------------------------------------------------------
1229
1230 end do
1231
1232 !--------------------------------------------------------------
1233
1234 stiff(1:nn*ndof, 1:nn*ndof) = tmpstiff(1:nn*ndof, 1:nn*ndof)
1235
1236 !--------------------------------------------------------------------
1237 end do
1238
1239 ! write(*,"(24E25.16)") stiff
1240
1241 !******************** Shell-Solid mixed analysis ********************
1242 ! mixglaf = 0 < natural shell (6 dof)
1243 ! mixglaf = 1 < mixed 361 (3*2 dof)*8 nod
1244 ! mixglaf = 2 < mixed 351 (3*2 dof)*6 nod
1245 if( mixflag == 1 )then
1246
1247 ! write(*,*) 'convert for shell-solid mixed analysis'
1248 sstable(1) = 1
1249 sstable(2) = 2
1250 sstable(3) = 3
1251 sstable(4) = 7
1252 sstable(5) = 8
1253 sstable(6) = 9
1254 sstable(7) = 13
1255 sstable(8) = 14
1256 sstable(9) = 15
1257 sstable(10)= 19
1258 sstable(11)= 20
1259 sstable(12)= 21
1260 sstable(13)= 4
1261 sstable(14)= 5
1262 sstable(15)= 6
1263 sstable(16)= 10
1264 sstable(17)= 11
1265 sstable(18)= 12
1266 sstable(19)= 16
1267 sstable(20)= 17
1268 sstable(21)= 18
1269 sstable(22)= 22
1270 sstable(23)= 23
1271 sstable(24)= 24
1272
1273 tmpstiff(1:nn*ndof, 1:nn*ndof) = stiff(1:nn*ndof, 1:nn*ndof)
1274
1275 do i = 1, nn*ndof
1276 do j = 1, nn*ndof
1277 stiff(i,j) = tmpstiff(sstable(i),sstable(j))
1278 enddo
1279 enddo
1280
1281 elseif( mixflag == 2 )then
1282 ! write(*,*) 'convert for shell-solid mixed analysis 351'
1283 sstable(1) = 1
1284 sstable(2) = 2
1285 sstable(3) = 3
1286 sstable(4) = 7
1287 sstable(5) = 8
1288 sstable(6) = 9
1289 sstable(7) = 13
1290 sstable(8) = 14
1291 sstable(9) = 15
1292 sstable(10)= 4
1293 sstable(11)= 5
1294 sstable(12)= 6
1295 sstable(13)= 10
1296 sstable(14)= 11
1297 sstable(15)= 12
1298 sstable(16)= 16
1299 sstable(17)= 17
1300 sstable(18)= 18
1301
1302 tmpstiff(1:nn*ndof, 1:nn*ndof) = stiff(1:nn*ndof, 1:nn*ndof)
1303
1304 do i = 1, nn*ndof
1305 do j = 1, nn*ndof
1306 stiff(i,j) = tmpstiff(sstable(i),sstable(j))
1307 enddo
1308 enddo
1309
1310 endif
1311
1312 return
1313
1314 !####################################################################
1315 end subroutine stf_shell_mitc
1316 !####################################################################
1317
1318
1319 !####################################################################
1321 (etype, nn, ndof, ecoord, gausses, edisp, &
1322 strain, stress, thick, zeta, n_layer, n_totlyr)
1323 !####################################################################
1324
1325 use mmechgauss
1327 use m_matmatrix
1328
1329 !--------------------------------------------------------------------
1330
1331 integer(kind = kint), intent(in) :: etype
1332 integer(kind = kint), intent(in) :: nn
1333 integer(kind = kint), intent(in) :: ndof
1334 real(kind = kreal), intent(in) :: ecoord(3, nn)
1335 type(tgaussstatus), intent(in) :: gausses(:)
1336 real(kind = kreal), intent(in) :: edisp(6, nn)
1337 real(kind = kreal), intent(out) :: strain(:,:)
1338 real(kind = kreal), intent(out) :: stress(:,:)
1339 real(kind = kreal), intent(in) :: thick
1340 real(kind = kreal), intent(in) :: zeta
1341
1342 !--------------------------------------------------------------------
1343
1344 integer :: i, m
1345 integer :: lx
1346 integer :: fetype
1347 integer :: ntying
1348 integer :: npoints_tying(3)
1349 integer :: it, ip
1350 integer :: na, nb
1351 integer :: n_layer, n_totlyr
1352
1353 real(kind = kreal) :: d(5, 5)
1354 real(kind = kreal) :: elem(3, nn)
1355 real(kind = kreal) :: xi_lx, eta_lx, zeta_ly
1356 real(kind = kreal) :: naturalcoord(2)
1357 real(kind = kreal) :: tpcoord(6, 2, 3)
1358 real(kind = kreal) :: nncoord(nn, 2)
1359 real(kind = kreal) :: shapefunc(nn)
1360 real(kind = kreal) :: shapederiv(nn, 2)
1361 real(kind = kreal) :: alpha
1362 real(kind = kreal) :: xxi_lx, eeta_lx
1363 real(kind = kreal) :: xxi_di(6, 3), eeta_di(6, 3)
1364 real(kind = kreal) :: h(nn, 3)
1365 real(kind = kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
1366 real(kind = kreal) :: v1_abs, v2_abs, v3_abs
1367 real(kind = kreal) :: a_over_2_v3(3, nn)
1368 real(kind = kreal) :: a_over_2_theta_cross_v3(3, nn)
1369 real(kind = kreal) :: u_rot(3, nn)
1370 real(kind = kreal) :: theta(3, nn)
1371 real(kind = kreal) :: dudxi(3), dudeta(3), dudzeta(3)
1372 real(kind = kreal) :: dudxi_rot(3, nn), dudeta_rot(3, nn), &
1373 dudzeta_rot(3, nn)
1374 real(kind = kreal) :: g1(3), g2(3), g3(3)
1375 real(kind = kreal) :: g3_abs
1376 real(kind = kreal) :: e_0(3)
1377 real(kind = kreal) :: cg1(3), cg2(3), cg3(3)
1378 real(kind = kreal) :: det
1379 real(kind = kreal) :: det_cg3(3)
1380 real(kind = kreal) :: det_inv
1381 real(kind = kreal) :: det_cg3_abs
1382 real(kind = kreal) :: e1_hat(3), e2_hat(3), e3_hat(3)
1383 real(kind = kreal) :: e1_hat_abs, e2_hat_abs
1384 real(kind = kreal) :: e11, e22, e12_2, e23_2, e31_2
1385 real(kind = kreal) :: e11_di(6, 3), e22_di(6, 3), &
1386 e12_di_2(6, 3), e23_di_2(6, 3), &
1387 e31_di_2(6, 3)
1388 real(kind = kreal) :: e(3, 3), ev(5)
1389 real(kind = kreal) :: s(3, 3), sv(5)
1390 real(kind = kreal) :: sumlyr
1391
1392
1393 zeta_ly = 0.0d0
1394
1395 !--------------------------------------------------------------------
1396
1397 ! for lamina stress
1398
1399 ! MITC4
1400 if( etype .EQ. fe_mitc4_shell ) then
1401
1402 fetype = fe_mitc4_shell
1403
1404 ntying = 1
1405 npoints_tying(1)= 4
1406
1407 ! MITC9
1408 else if( etype .EQ. fe_mitc9_shell ) then
1409
1410 fetype = fe_mitc9_shell
1411
1412 ntying = 3
1413 npoints_tying(1)= 6
1414 npoints_tying(2)= 6
1415 npoints_tying(3)= 4
1416
1417 ! MITC3
1418 else if( etype .EQ. fe_mitc3_shell ) then
1419
1420 fetype = fe_mitc3_shell
1421
1422 ntying = 1
1423 npoints_tying(1)= 3
1424
1425 end if
1426
1427 !--------------------------------------------------------------------
1428
1429 elem(:, :) = ecoord(:, :)
1430
1431 !--------------------------------------------------------------------
1432
1433 do na = 1, nn
1434
1435 theta(1, na) = edisp(4, na)
1436 theta(2, na) = edisp(5, na)
1437 theta(3, na) = edisp(6, na)
1438
1439 end do
1440
1441 !-------------------------------------------------------------------
1442
1443 ! xi-coordinate at a node in a local element
1444 ! eta-coordinate at a node in a local element
1445 call getnodalnaturalcoord(fetype, nncoord)
1446
1447 !-------------------------------------------------------------------
1448
1449 ! MITC4
1450 if( etype .EQ. fe_mitc4_shell ) then
1451
1452 !--------------------------------------------------------
1453
1454 ! xi-coordinate at a tying pont in a local element
1455 tpcoord(1, 1, 1) = 0.0d0
1456 tpcoord(2, 1, 1) = 1.0d0
1457 tpcoord(3, 1, 1) = 0.0d0
1458 tpcoord(4, 1, 1) = -1.0d0
1459 ! eta-coordinate at a tying point in a local element
1460 tpcoord(1, 2, 1) = -1.0d0
1461 tpcoord(2, 2, 1) = 0.0d0
1462 tpcoord(3, 2, 1) = 1.0d0
1463 tpcoord(4, 2, 1) = 0.0d0
1464
1465 !--------------------------------------------------------
1466
1467 ! MITC9
1468 else if( etype .EQ. fe_mitc9_shell ) then
1469
1470 !--------------------------------------------------------
1471
1472 ! xi-coordinate at a tying point in a local element
1473 tpcoord(1, 1, 1) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
1474 tpcoord(2, 1, 1) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
1475 tpcoord(3, 1, 1) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
1476 tpcoord(4, 1, 1) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
1477 tpcoord(5, 1, 1) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
1478 tpcoord(6, 1, 1) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
1479 ! eta-coordinate at a tying point in a local element
1480 tpcoord(1, 2, 1) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
1481 tpcoord(2, 2, 1) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
1482 tpcoord(3, 2, 1) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
1483 tpcoord(4, 2, 1) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
1484 tpcoord(5, 2, 1) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
1485 tpcoord(6, 2, 1) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
1486
1487 ! xi-coordinate at a tying point in a local element
1488 tpcoord(1, 1, 2) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
1489 tpcoord(2, 1, 2) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
1490 tpcoord(3, 1, 2) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
1491 tpcoord(4, 1, 2) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
1492 tpcoord(5, 1, 2) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
1493 tpcoord(6, 1, 2) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
1494 ! eta-coordinate at a tying point in a local element
1495 tpcoord(1, 2, 2) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
1496 tpcoord(2, 2, 2) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
1497 tpcoord(3, 2, 2) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
1498 tpcoord(4, 2, 2) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
1499 tpcoord(5, 2, 2) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
1500 tpcoord(6, 2, 2) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
1501
1502 ! xi-coordinate at a tying point in a local element
1503 tpcoord(1, 1, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
1504 tpcoord(2, 1, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
1505 tpcoord(3, 1, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
1506 tpcoord(4, 1, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
1507 ! eta-coordinate at a tying point in a local element
1508 tpcoord(1, 2, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
1509 tpcoord(2, 2, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
1510 tpcoord(3, 2, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
1511 tpcoord(4, 2, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
1512
1513 !--------------------------------------------------------
1514
1515 ! Xi-coordinate at a tying point in a local element
1516 xxi_di(1, 1) = -1.0d0
1517 xxi_di(2, 1) = 1.0d0
1518 xxi_di(3, 1) = 1.0d0
1519 xxi_di(4, 1) = -1.0d0
1520 xxi_di(5, 1) = 1.0d0
1521 xxi_di(6, 1) = -1.0d0
1522 ! Eta-coordinate at a tying point in a local element
1523 eeta_di(1, 1) = -1.0d0
1524 eeta_di(2, 1) = -1.0d0
1525 eeta_di(3, 1) = 1.0d0
1526 eeta_di(4, 1) = 1.0d0
1527 eeta_di(5, 1) = 0.0d0
1528 eeta_di(6, 1) = 0.0d0
1529
1530 ! Xi-coordinate at a tying point in a local element
1531 xxi_di(1, 2) = -1.0d0
1532 xxi_di(2, 2) = 0.0d0
1533 xxi_di(3, 2) = 1.0d0
1534 xxi_di(4, 2) = 1.0d0
1535 xxi_di(5, 2) = 0.0d0
1536 xxi_di(6, 2) = -1.0d0
1537 ! Eta-coordinate at a tying point in a local element
1538 eeta_di(1, 2) = -1.0d0
1539 eeta_di(2, 2) = -1.0d0
1540 eeta_di(3, 2) = -1.0d0
1541 eeta_di(4, 2) = 1.0d0
1542 eeta_di(5, 2) = 1.0d0
1543 eeta_di(6, 2) = 1.0d0
1544
1545 !--------------------------------------------------------
1546
1547 ! MITC3
1548 else if( etype .EQ. fe_mitc3_shell ) then
1549
1550 !--------------------------------------------------------
1551
1552 ! xi-coordinate at a tying point in a local element
1553 tpcoord(1, 1, 1) = 0.5d0
1554 tpcoord(2, 1, 1) = 0.0d0
1555 tpcoord(3, 1, 1) = 0.5d0
1556 ! eta-coordinate at a tying point in a local element
1557 tpcoord(1, 2, 1) = 0.0d0
1558 tpcoord(2, 2, 1) = 0.5d0
1559 tpcoord(3, 2, 1) = 0.5d0
1560
1561 !--------------------------------------------------------
1562
1563 end if
1564
1565 !--------------------------------------------------------------------
1566
1567 ! xi-coordinate at the center point in a local element
1568 ! eta-coordinate at the center point in a local element
1569 naturalcoord(1) = 0.0d0
1570 naturalcoord(2) = 0.0d0
1571
1572 call getshapederiv(fetype, naturalcoord, shapederiv)
1573
1574 !--------------------------------------------------------------
1575
1576 ! Covariant basis vector
1577 do i = 1, 3
1578
1579 g1(i) = 0.0d0
1580
1581 do na = 1, nn
1582
1583 g1(i) = g1(i)+shapederiv(na, 1) &
1584 *elem(i, na)
1585
1586 end do
1587
1588 end do
1589
1590 e_0(1) = g1(1)
1591 e_0(2) = g1(2)
1592 e_0(3) = g1(3)
1593
1594 !--------------------------------------------------------------
1595
1596 do nb = 1, nn
1597
1598 !--------------------------------------------------------
1599
1600 naturalcoord(1) = nncoord(nb, 1)
1601 naturalcoord(2) = nncoord(nb, 2)
1602
1603 call getshapederiv(fetype, naturalcoord, shapederiv)
1604
1605 !--------------------------------------------------------
1606
1607 ! Covariant basis vector
1608 do i = 1, 3
1609
1610 g1(i) = 0.0d0
1611 g2(i) = 0.0d0
1612
1613 do na = 1, nn
1614
1615 g1(i) = g1(i)+shapederiv(na, 1) &
1616 *elem(i, na)
1617 g2(i) = g2(i)+shapederiv(na, 2) &
1618 *elem(i, na)
1619
1620 end do
1621
1622 end do
1623
1624 !--------------------------------------------------------
1625
1626 det_cg3(1) = g1(2)*g2(3)-g1(3)*g2(2)
1627 det_cg3(2) = g1(3)*g2(1)-g1(1)*g2(3)
1628 det_cg3(3) = g1(1)*g2(2)-g1(2)*g2(1)
1629
1630 det_cg3_abs = dsqrt( det_cg3(1)*det_cg3(1) &
1631 +det_cg3(2)*det_cg3(2) &
1632 +det_cg3(3)*det_cg3(3) )
1633
1634 v3(1, nb) = det_cg3(1)/det_cg3_abs
1635 v3(2, nb) = det_cg3(2)/det_cg3_abs
1636 v3(3, nb) = det_cg3(3)/det_cg3_abs
1637
1638 !--------------------------------------------------------
1639
1640 v2(1, nb) = v3(2, nb)*e_0(3)-v3(3, nb)*e_0(2)
1641 v2(2, nb) = v3(3, nb)*e_0(1)-v3(1, nb)*e_0(3)
1642 v2(3, nb) = v3(1, nb)*e_0(2)-v3(2, nb)*e_0(1)
1643
1644 v2_abs = dsqrt( v2(1, nb)*v2(1, nb) &
1645 +v2(2, nb)*v2(2, nb) &
1646 +v2(3, nb)*v2(3, nb) )
1647
1648 if( v2_abs .GT. 1.0d-15 ) then
1649
1650 v2(1, nb) = v2(1, nb)/v2_abs
1651 v2(2, nb) = v2(2, nb)/v2_abs
1652 v2(3, nb) = v2(3, nb)/v2_abs
1653
1654 v1(1, nb) = v2(2, nb)*v3(3, nb) &
1655 -v2(3, nb)*v3(2, nb)
1656 v1(2, nb) = v2(3, nb)*v3(1, nb) &
1657 -v2(1, nb)*v3(3, nb)
1658 v1(3, nb) = v2(1, nb)*v3(2, nb) &
1659 -v2(2, nb)*v3(1, nb)
1660
1661 v1_abs = dsqrt( v1(1, nb)*v1(1, nb) &
1662 +v1(2, nb)*v1(2, nb) &
1663 +v1(3, nb)*v1(3, nb) )
1664
1665 v1(1, nb) = v1(1, nb)/v1_abs
1666 v1(2, nb) = v1(2, nb)/v1_abs
1667 v1(3, nb) = v1(3, nb)/v1_abs
1668
1669 else
1670
1671 v1(1, nb) = 0.0d0
1672 v1(2, nb) = 0.0d0
1673 v1(3, nb) = -1.0d0
1674
1675 v2(1, nb) = 0.0d0
1676 v2(2, nb) = 1.0d0
1677 v2(3, nb) = 0.0d0
1678
1679 end if
1680
1681 !--------------------------------------------------------
1682
1683 v3(1, nb) = v1(2, nb)*v2(3, nb) &
1684 -v1(3, nb)*v2(2, nb)
1685 v3(2, nb) = v1(3, nb)*v2(1, nb) &
1686 -v1(1, nb)*v2(3, nb)
1687 v3(3, nb) = v1(1, nb)*v2(2, nb) &
1688 -v1(2, nb)*v2(1, nb)
1689
1690 v3_abs = dsqrt( v3(1, nb)*v3(1, nb) &
1691 +v3(2, nb)*v3(2, nb) &
1692 +v3(3, nb)*v3(3, nb) )
1693
1694 v3(1, nb) = v3(1, nb)/v3_abs
1695 v3(2, nb) = v3(2, nb)/v3_abs
1696 v3(3, nb) = v3(3, nb)/v3_abs
1697
1698 !--------------------------------------------------------
1699
1700 a_over_2_v3(1, nb) = 0.5d0*thick*v3(1, nb)
1701 a_over_2_v3(2, nb) = 0.5d0*thick*v3(2, nb)
1702 a_over_2_v3(3, nb) = 0.5d0*thick*v3(3, nb)
1703
1704 !--------------------------------------------------------
1705
1706 a_over_2_theta_cross_v3(1, nb) &
1707 = theta(2, nb)*a_over_2_v3(3, nb) &
1708 -theta(3, nb)*a_over_2_v3(2, nb)
1709 a_over_2_theta_cross_v3(2, nb) &
1710 = theta(3, nb)*a_over_2_v3(1, nb) &
1711 -theta(1, nb)*a_over_2_v3(3, nb)
1712 a_over_2_theta_cross_v3(3, nb) &
1713 = theta(1, nb)*a_over_2_v3(2, nb) &
1714 -theta(2, nb)*a_over_2_v3(1, nb)
1715
1716 !--------------------------------------------------------
1717
1718 end do
1719
1720 !--------------------------------------------------------------------
1721 ! Modified stress in laminated shell
1722 !--------------------------------------------------------------------
1723 ! MITC4
1724 if( etype .EQ. fe_mitc4_shell ) then
1725
1726 zeta_ly = 0.0d0
1727
1728 ! MITC9
1729 else if( etype .EQ. fe_mitc9_shell ) then
1730
1731 zeta_ly = zeta
1732
1733 ! MITC3
1734 else if( etype .EQ. fe_mitc3_shell ) then
1735
1736 zeta_ly = 0.0d0
1737
1738 end if
1739
1740 !---------------------------------------------------------
1741
1742 do it = 1, ntying
1743
1744 do ip = 1, npoints_tying(it)
1745
1746 !-------------------------------------------------
1747
1748 naturalcoord(1) = tpcoord(ip, 1, it)
1749 naturalcoord(2) = tpcoord(ip, 2, it)
1750
1751 call getshapefunc(fetype, naturalcoord, shapefunc)
1752
1753 call getshapederiv(fetype, naturalcoord, shapederiv)
1754
1755 !-------------------------------------------------
1756
1757 do na = 1, nn
1758
1759 do i = 1, 3
1760
1761 u_rot(i, na) &
1762 = shapefunc(na) &
1763 *( zeta_ly*a_over_2_v3(i, na) )
1764
1765 dudxi_rot(i, na) &
1766 = shapederiv(na, 1) &
1767 *( zeta_ly*a_over_2_v3(i, na) )
1768 dudeta_rot(i, na) &
1769 = shapederiv(na, 2) &
1770 *( zeta_ly*a_over_2_v3(i, na) )
1771 dudzeta_rot(i, na) &
1772 = shapefunc(na) &
1773 *( a_over_2_v3(i, na) )
1774
1775 end do
1776
1777 end do
1778
1779 !-------------------------------------------------
1780
1781 ! Covariant basis vector
1782 do i = 1, 3
1783
1784 g1(i) = 0.0d0
1785 g2(i) = 0.0d0
1786 g3(i) = 0.0d0
1787
1788 do na = 1, nn
1789
1790 g1(i) = g1(i)+shapederiv(na, 1) &
1791 *elem(i, na) &
1792 +dudxi_rot(i, na)
1793 g2(i) = g2(i)+shapederiv(na, 2) &
1794 *elem(i, na) &
1795 +dudeta_rot(i, na)
1796 g3(i) = g3(i)+dudzeta_rot(i, na)
1797
1798 end do
1799
1800 end do
1801
1802 !---------------------------------------------
1803
1804 do i = 1, 3
1805
1806 dudxi(i) = 0.0d0
1807 dudeta(i) = 0.0d0
1808 dudzeta(i) = 0.0d0
1809
1810 do na = 1, nn
1811
1812 dudxi(i) &
1813 = dudxi(i) &
1814 +shapederiv(na, 1) &
1815 *( edisp(i, na) &
1816 +zeta_ly*a_over_2_theta_cross_v3(i, na) )
1817 dudeta(i) &
1818 = dudeta(i) &
1819 +shapederiv(na, 2) &
1820 *( edisp(i, na) &
1821 +zeta_ly*a_over_2_theta_cross_v3(i, na) )
1822 dudzeta(i) &
1823 = dudzeta(i) &
1824 +shapefunc(na) &
1825 *( a_over_2_theta_cross_v3(i, na) )
1826
1827 end do
1828
1829 end do
1830
1831 !---------------------------------------------
1832
1833 ! Infinitesimal strain tensor
1834 e11_di(ip, it) &
1835 = 0.5d0 &
1836 *( ( g1(1)*dudxi(1) +dudxi(1) *g1(1) ) &
1837 +( g1(2)*dudxi(2) +dudxi(2) *g1(2) ) &
1838 +( g1(3)*dudxi(3) +dudxi(3) *g1(3) ) )
1839 e22_di(ip, it) &
1840 = 0.5d0 &
1841 *( ( g2(1)*dudeta(1)+dudeta(1)*g2(1) ) &
1842 +( g2(2)*dudeta(2)+dudeta(2)*g2(2) ) &
1843 +( g2(3)*dudeta(3)+dudeta(3)*g2(3) ) )
1844 e12_2 &
1845 = ( g1(1)*dudeta(1) +dudxi(1) *g2(1) ) &
1846 +( g1(2)*dudeta(2) +dudxi(2) *g2(2) ) &
1847 +( g1(3)*dudeta(3) +dudxi(3) *g2(3) )
1848 e23_di_2(ip, it) &
1849 = ( g2(1)*dudzeta(1)+dudeta(1) *g3(1) ) &
1850 +( g2(2)*dudzeta(2)+dudeta(2) *g3(2) ) &
1851 +( g2(3)*dudzeta(3)+dudeta(3) *g3(3) )
1852 e31_di_2(ip, it) &
1853 = ( g3(1)*dudxi(1) +dudzeta(1)*g1(1) ) &
1854 +( g3(2)*dudxi(2) +dudzeta(2)*g1(2) ) &
1855 +( g3(3)*dudxi(3) +dudzeta(3)*g1(3) )
1856
1857 !-------------------------------------------------
1858
1859 end do
1860
1861 end do
1862
1863 !--------------------------------------------------------
1864
1865 sumlyr = 0.0d0
1866 do m = 1, n_layer
1867 sumlyr = sumlyr +2*gausses(1)%pMaterial%shell_var(m)%weight
1868 end do
1869 zeta_ly = -1 +sumlyr -gausses(1)%pMaterial%shell_var(n_layer)%weight*(1-zeta)
1870
1871 !--------------------------------------------------------
1872
1873 do lx = 1, nn
1874
1875 !--------------------------------------------------
1876
1877 naturalcoord(1) = nncoord(lx, 1)
1878 naturalcoord(2) = nncoord(lx, 2)
1879
1880 xi_lx = naturalcoord(1)
1881 eta_lx = naturalcoord(2)
1882
1883 call getshapefunc(fetype, naturalcoord, shapefunc)
1884
1885 call getshapederiv(fetype, naturalcoord, shapederiv)
1886
1887 !--------------------------------------------------
1888
1889 do na = 1, nn
1890
1891 do i = 1, 3
1892
1893 u_rot(i, na) &
1894 = shapefunc(na) &
1895 *( zeta_ly*a_over_2_v3(i, na) )
1896
1897 dudxi_rot(i, na) &
1898 = shapederiv(na, 1) &
1899 *( zeta_ly*a_over_2_v3(i, na) )
1900 dudeta_rot(i, na) &
1901 = shapederiv(na, 2) &
1902 *( zeta_ly*a_over_2_v3(i, na) )
1903 dudzeta_rot(i, na) &
1904 = shapefunc(na) &
1905 *( a_over_2_v3(i, na) )
1906
1907 end do
1908
1909 end do
1910
1911 !--------------------------------------------------
1912
1913 ! Covariant basis vector
1914 do i = 1, 3
1915
1916 g1(i) = 0.0d0
1917 g2(i) = 0.0d0
1918 g3(i) = 0.0d0
1919
1920 do na = 1, nn
1921
1922 g1(i) = g1(i)+shapederiv(na, 1) &
1923 *elem(i, na) &
1924 +dudxi_rot(i, na)
1925 g2(i) = g2(i)+shapederiv(na, 2) &
1926 *elem(i, na) &
1927 +dudeta_rot(i, na)
1928 g3(i) = g3(i)+dudzeta_rot(i, na)
1929
1930 end do
1931 end do
1932
1933 !--------------------------------------------------
1934
1935 ! Jacobian
1936 det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
1937 +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
1938 +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )
1939
1940 if(det == 0.0d0) then
1941 write(*,*)"ERROR:LIB Shell in l2009 Not Jacobian"
1942 stop
1943 endif
1944
1945 det_inv = 1.0d0/det
1946
1947 !--------------------------------------------------
1948
1949 ! Contravariant basis vector
1950 cg1(1) = det_inv &
1951 *( g2(2)*g3(3)-g2(3)*g3(2) )
1952 cg1(2) = det_inv &
1953 *( g2(3)*g3(1)-g2(1)*g3(3) )
1954 cg1(3) = det_inv &
1955 *( g2(1)*g3(2)-g2(2)*g3(1) )
1956 cg2(1) = det_inv &
1957 *( g3(2)*g1(3)-g3(3)*g1(2) )
1958 cg2(2) = det_inv &
1959 *( g3(3)*g1(1)-g3(1)*g1(3) )
1960 cg2(3) = det_inv &
1961 *( g3(1)*g1(2)-g3(2)*g1(1) )
1962 cg3(1) = det_inv &
1963 *( g1(2)*g2(3)-g1(3)*g2(2) )
1964 cg3(2) = det_inv &
1965 *( g1(3)*g2(1)-g1(1)*g2(3) )
1966 cg3(3) = det_inv &
1967 *( g1(1)*g2(2)-g1(2)*g2(1) )
1968
1969 !--------------------------------------------------
1970
1971 g3_abs = dsqrt( g3(1)*g3(1) &
1972 +g3(2)*g3(2) &
1973 +g3(3)*g3(3) )
1974
1975 !--------------------------------------------------
1976
1977 ! Orthonormal vectors
1978
1979 e3_hat(1) = g3(1)/g3_abs
1980 e3_hat(2) = g3(2)/g3_abs
1981 e3_hat(3) = g3(3)/g3_abs
1982
1983 e1_hat(1) = g2(2)*e3_hat(3) &
1984 -g2(3)*e3_hat(2)
1985 e1_hat(2) = g2(3)*e3_hat(1) &
1986 -g2(1)*e3_hat(3)
1987 e1_hat(3) = g2(1)*e3_hat(2) &
1988 -g2(2)*e3_hat(1)
1989 e1_hat_abs = dsqrt( e1_hat(1)*e1_hat(1) &
1990 +e1_hat(2)*e1_hat(2) &
1991 +e1_hat(3)*e1_hat(3) )
1992 e1_hat(1) = e1_hat(1)/e1_hat_abs
1993 e1_hat(2) = e1_hat(2)/e1_hat_abs
1994 e1_hat(3) = e1_hat(3)/e1_hat_abs
1995
1996 e2_hat(1) = e3_hat(2)*e1_hat(3) &
1997 -e3_hat(3)*e1_hat(2)
1998 e2_hat(2) = e3_hat(3)*e1_hat(1) &
1999 -e3_hat(1)*e1_hat(3)
2000 e2_hat(3) = e3_hat(1)*e1_hat(2) &
2001 -e3_hat(2)*e1_hat(1)
2002 e2_hat_abs = dsqrt( e2_hat(1)*e2_hat(1) &
2003 +e2_hat(2)*e2_hat(2) &
2004 +e2_hat(3)*e2_hat(3) )
2005 e2_hat(1) = e2_hat(1)/e2_hat_abs
2006 e2_hat(2) = e2_hat(2)/e2_hat_abs
2007 e2_hat(3) = e2_hat(3)/e2_hat_abs
2008
2009 !--------------------------------------------------
2010
2011 do i = 1, 3
2012
2013 dudxi(i) = 0.0d0
2014 dudeta(i) = 0.0d0
2015 dudzeta(i) = 0.0d0
2016
2017 do na = 1, nn
2018
2019 dudxi(i) &
2020 = dudxi(i) &
2021 +shapederiv(na, 1) &
2022 *( edisp(i, na) &
2023 +zeta_ly*a_over_2_theta_cross_v3(i, na) )
2024 dudeta(i) &
2025 = dudeta(i) &
2026 +shapederiv(na, 2) &
2027 *( edisp(i, na) &
2028 +zeta_ly*a_over_2_theta_cross_v3(i, na) )
2029 dudzeta(i) &
2030 = dudzeta(i) &
2031 +shapefunc(na) &
2032 *( a_over_2_theta_cross_v3(i, na) )
2033
2034 end do
2035
2036 end do
2037
2038 !--------------------------------------------------
2039
2040 ! Infinitesimal strain tensor
2041 e11 &
2042 = 0.5d0 &
2043 *( ( g1(1)*dudxi(1) +dudxi(1) *g1(1) ) &
2044 +( g1(2)*dudxi(2) +dudxi(2) *g1(2) ) &
2045 +( g1(3)*dudxi(3) +dudxi(3) *g1(3) ) )
2046 e22 &
2047 = 0.5d0 &
2048 *( ( g2(1)*dudeta(1)+dudeta(1)*g2(1) ) &
2049 +( g2(2)*dudeta(2)+dudeta(2)*g2(2) ) &
2050 +( g2(3)*dudeta(3)+dudeta(3)*g2(3) ) )
2051 e12_2 &
2052 = ( g1(1)*dudeta(1) +dudxi(1) *g2(1) ) &
2053 +( g1(2)*dudeta(2) +dudxi(2) *g2(2) ) &
2054 +( g1(3)*dudeta(3) +dudxi(3) *g2(3) )
2055 e23_2 &
2056 = ( g2(1)*dudzeta(1)+dudeta(1) *g3(1) ) &
2057 +( g2(2)*dudzeta(2)+dudeta(2) *g3(2) ) &
2058 +( g2(3)*dudzeta(3)+dudeta(3) *g3(3) )
2059 e31_2 &
2060 = ( g3(1)*dudxi(1) +dudzeta(1)*g1(1) ) &
2061 +( g3(2)*dudxi(2) +dudzeta(2)*g1(2) ) &
2062 +( g3(3)*dudxi(3) +dudzeta(3)*g1(3) )
2063
2064 ! MITC4
2065 if( etype .EQ. fe_mitc4_shell ) then
2066
2067 e23_2 = 0.0d0
2068 e31_2 = 0.0d0
2069
2070 ! e23_as_2
2071 e23_2 &
2072 = 0.5d0*( 1.0d0-xi_lx )*e23_di_2(4, 1) &
2073 +0.5d0*( 1.0d0+xi_lx )*e23_di_2(2, 1)
2074 ! e31_as_2
2075 e31_2 &
2076 = 0.5d0*( 1.0d0-eta_lx )*e31_di_2(1, 1) &
2077 +0.5d0*( 1.0d0+eta_lx )*e31_di_2(3, 1)
2078
2079 ! MITC9
2080 else if( etype .EQ. fe_mitc9_shell ) then
2081
2082 xxi_lx = xi_lx /dsqrt( 1.0d0/3.0d0 )
2083 eeta_lx = eta_lx/dsqrt( 3.0d0/5.0d0 )
2084
2085 do ip = 1, npoints_tying(1)
2086
2087 h(ip, 1) &
2088 = ( 0.5d0*( 1.0d0+xxi_di(ip, 1)*xxi_lx ) ) &
2089 *( ( 0.5d0*eeta_di(ip, 1)*eeta_lx ) &
2090 *( 1.0d0+eeta_di(ip, 1)*eeta_lx ) &
2091 +( 1.0d0-eeta_di(ip, 1)*eeta_di(ip, 1) ) &
2092 *( 1.0d0-eeta_lx*eeta_lx ) )
2093
2094 end do
2095
2096 xxi_lx = xi_lx /dsqrt( 3.0d0/5.0d0 )
2097 eeta_lx = eta_lx/dsqrt( 1.0d0/3.0d0 )
2098
2099 do ip = 1, npoints_tying(2)
2100
2101 h(ip, 2) &
2102 = ( ( 0.5d0*xxi_di(ip, 2) *xxi_lx ) &
2103 *( 1.0d0+xxi_di(ip, 2) *xxi_lx ) &
2104 +( 1.0d0-xxi_di(ip, 2) *xxi_di(ip, 2) ) &
2105 *( 1.0d0-xxi_lx*xxi_lx ) ) &
2106 *( 0.5d0*( 1.0d0+eeta_di(ip, 2)*eeta_lx ) )
2107
2108 end do
2109
2110 xxi_lx = xi_lx /dsqrt( 1.0d0/3.0d0 )
2111 eeta_lx = eta_lx/dsqrt( 1.0d0/3.0d0 )
2112
2113 do ip = 1, npoints_tying(3)
2114
2115 h(ip, 3) &
2116 = ( 0.5d0*( 1.0d0+xxi_di(ip, 1)*xxi_lx ) ) &
2117 *( 0.5d0*( 1.0d0+eeta_di(ip, 1)*eeta_lx ) )
2118
2119 end do
2120
2121 e11 = 0.0d0
2122 e31_2 = 0.0d0
2123
2124 ! e11_as, e31_as_2
2125 do ip = 1, npoints_tying(1)
2126
2127 e11 = e11 +h(ip, 1)*e11_di(ip, 1)
2128 e31_2 = e31_2+h(ip, 1)*e31_di_2(ip, 1)
2129
2130 end do
2131
2132 e22 = 0.0d0
2133 e23_2 = 0.0d0
2134
2135 ! e22_as, e23_as_2
2136 do ip = 1, npoints_tying(2)
2137
2138 e22 = e22 +h(ip, 2)*e22_di(ip, 2)
2139 e23_2 = e23_2+h(ip, 2)*e23_di_2(ip, 2)
2140
2141 end do
2142
2143 e12_2 = 0.0d0
2144
2145 ! e12_as_2
2146 do ip = 1, npoints_tying(3)
2147
2148 e12_2 = e12_2+h(ip, 3)*e12_di_2(ip, 3)
2149
2150 end do
2151
2152 ! MITC3
2153 else if( etype .EQ. fe_mitc3_shell ) then
2154
2155 e23_2 = 0.0d0
2156 e31_2 = 0.0d0
2157
2158 ! e23_as_2
2159 e23_2 &
2160 = ( 1.0d0-xi_lx )*e23_di_2(2, 1) &
2161 +xi_lx *e31_di_2(1, 1) &
2162 +xi_lx *( e23_di_2(3, 1)-e31_di_2(3, 1) )
2163 ! e31_as_2
2164 e31_2 &
2165 = eta_lx*e23_di_2(2, 1) &
2166 +( 1.0d0-eta_lx )*e31_di_2(1, 1) &
2167 -eta_lx*( e23_di_2(3, 1)-e31_di_2(3, 1) )
2168
2169 end if
2170
2171 !--------------------------------------------------
2172
2173 ! { E } vector
2174 ev(1) = e11
2175 ev(2) = e22
2176 ev(3) = e12_2
2177 ev(4) = e23_2
2178 ev(5) = e31_2
2179
2180 ! Infinitesimal strain tensor
2181 ! [ E ] matrix
2182 e(1, 1) = ev(1)
2183 e(2, 2) = ev(2)
2184 e(3, 3) = 0.0d0
2185 e(1, 2) = 0.5d0*ev(3)
2186 e(2, 1) = 0.5d0*ev(3)
2187 e(2, 3) = 0.5d0*ev(4)
2188 e(3, 2) = 0.5d0*ev(4)
2189 e(3, 1) = 0.5d0*ev(5)
2190 e(1, 3) = 0.5d0*ev(5)
2191
2192 !--------------------------------------------------
2193 ! write(*,*) 'Stress_n_layer', n_layer
2194 call matlmatrix_shell &
2195 (gausses(lx), shell, d, &
2196 e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
2197 alpha, n_layer)
2198
2199 !--------------------------------------------------
2200
2201 sv = matmul( d, ev )
2202
2203 ! Infinitesimal stress tensor
2204 ! [ S ] matrix
2205 s(1, 1) = sv(1)
2206 s(2, 2) = sv(2)
2207 s(3, 3) = 0.0d0
2208 s(1, 2) = sv(3)
2209 s(2, 1) = sv(3)
2210 s(2, 3) = sv(4)
2211 s(3, 2) = sv(4)
2212 s(3, 1) = sv(5)
2213 s(1, 3) = sv(5)
2214
2215 !------------------------------------------------
2216
2217 stress(lx,1) &
2218 = ( s(1, 1)*g1(1)*g1(1) &
2219 +s(1, 2)*g1(1)*g2(1) &
2220 +s(1, 3)*g1(1)*g3(1) &
2221 +s(2, 1)*g2(1)*g1(1) &
2222 +s(2, 2)*g2(1)*g2(1) &
2223 +s(2, 3)*g2(1)*g3(1) &
2224 +s(3, 1)*g3(1)*g1(1) &
2225 +s(3, 2)*g3(1)*g2(1) )
2226 stress(lx,2) &
2227 = ( s(1, 1)*g1(2)*g1(2) &
2228 +s(1, 2)*g1(2)*g2(2) &
2229 +s(1, 3)*g1(2)*g3(2) &
2230 +s(2, 1)*g2(2)*g1(2) &
2231 +s(2, 2)*g2(2)*g2(2) &
2232 +s(2, 3)*g2(2)*g3(2) &
2233 +s(3, 1)*g3(2)*g1(2) &
2234 +s(3, 2)*g3(2)*g2(2) )
2235 stress(lx,3) &
2236 = ( s(1, 1)*g1(3)*g1(3) &
2237 +s(1, 2)*g1(3)*g2(3) &
2238 +s(1, 3)*g1(3)*g3(3) &
2239 +s(2, 1)*g2(3)*g1(3) &
2240 +s(2, 2)*g2(3)*g2(3) &
2241 +s(2, 3)*g2(3)*g3(3) &
2242 +s(3, 1)*g3(3)*g1(3) &
2243 +s(3, 2)*g3(3)*g2(3) )
2244 stress(lx,4) &
2245 = ( s(1, 1)*g1(1)*g1(2) &
2246 +s(1, 2)*g1(1)*g2(2) &
2247 +s(1, 3)*g1(1)*g3(2) &
2248 +s(2, 1)*g2(1)*g1(2) &
2249 +s(2, 2)*g2(1)*g2(2) &
2250 +s(2, 3)*g2(1)*g3(2) &
2251 +s(3, 1)*g3(1)*g1(2) &
2252 +s(3, 2)*g3(1)*g2(2) )
2253 stress(lx,5) &
2254 = ( s(1, 1)*g1(2)*g1(3) &
2255 +s(1, 2)*g1(2)*g2(3) &
2256 +s(1, 3)*g1(2)*g3(3) &
2257 +s(2, 1)*g2(2)*g1(3) &
2258 +s(2, 2)*g2(2)*g2(3) &
2259 +s(2, 3)*g2(2)*g3(3) &
2260 +s(3, 1)*g3(2)*g1(3) &
2261 +s(3, 2)*g3(2)*g2(3) )
2262 stress(lx,6) &
2263 = ( s(1, 1)*g1(3)*g1(1) &
2264 +s(1, 2)*g1(3)*g2(1) &
2265 +s(1, 3)*g1(3)*g3(1) &
2266 +s(2, 1)*g2(3)*g1(1) &
2267 +s(2, 2)*g2(3)*g2(1) &
2268 +s(2, 3)*g2(3)*g3(1) &
2269 +s(3, 1)*g3(3)*g1(1) &
2270 +s(3, 2)*g3(3)*g2(1) )
2271
2272 strain(lx,1) &
2273 = ( e(1, 1)*cg1(1)*cg1(1) &
2274 +e(1, 2)*cg1(1)*cg2(1) &
2275 +e(1, 3)*cg1(1)*cg3(1) &
2276 +e(2, 1)*cg2(1)*cg1(1) &
2277 +e(2, 2)*cg2(1)*cg2(1) &
2278 +e(2, 3)*cg2(1)*cg3(1) &
2279 +e(3, 1)*cg3(1)*cg1(1) &
2280 +e(3, 2)*cg3(1)*cg2(1) )
2281 strain(lx,2) &
2282 = ( e(1, 1)*cg1(2)*cg1(2) &
2283 +e(1, 2)*cg1(2)*cg2(2) &
2284 +e(1, 3)*cg1(2)*cg3(2) &
2285 +e(2, 1)*cg2(2)*cg1(2) &
2286 +e(2, 2)*cg2(2)*cg2(2) &
2287 +e(2, 3)*cg2(2)*cg3(2) &
2288 +e(3, 1)*cg3(2)*cg1(2) &
2289 +e(3, 2)*cg3(2)*cg2(2) )
2290 strain(lx,3) &
2291 = ( e(1, 1)*cg1(3)*cg1(3) &
2292 +e(1, 2)*cg1(3)*cg2(3) &
2293 +e(1, 3)*cg1(3)*cg3(3) &
2294 +e(2, 1)*cg2(3)*cg1(3) &
2295 +e(2, 2)*cg2(3)*cg2(3) &
2296 +e(2, 3)*cg2(3)*cg3(3) &
2297 +e(3, 1)*cg3(3)*cg1(3) &
2298 +e(3, 2)*cg3(3)*cg2(3) )
2299 strain(lx,4) &
2300 = ( e(1, 1)*cg1(1)*cg1(2) &
2301 +e(1, 2)*cg1(1)*cg2(2) &
2302 +e(1, 3)*cg1(1)*cg3(2) &
2303 +e(2, 1)*cg2(1)*cg1(2) &
2304 +e(2, 2)*cg2(1)*cg2(2) &
2305 +e(2, 3)*cg2(1)*cg3(2) &
2306 +e(3, 1)*cg3(1)*cg1(2) &
2307 +e(3, 2)*cg3(1)*cg2(2) )
2308 strain(lx,5) &
2309 = ( e(1, 1)*cg1(2)*cg1(3) &
2310 +e(1, 2)*cg1(2)*cg2(3) &
2311 +e(1, 3)*cg1(2)*cg3(3) &
2312 +e(2, 1)*cg2(2)*cg1(3) &
2313 +e(2, 2)*cg2(2)*cg2(3) &
2314 +e(2, 3)*cg2(2)*cg3(3) &
2315 +e(3, 1)*cg3(2)*cg1(3) &
2316 +e(3, 2)*cg3(2)*cg2(3) )
2317 strain(lx,6) &
2318 = ( e(1, 1)*cg1(3)*cg1(1) &
2319 +e(1, 2)*cg1(3)*cg2(1) &
2320 +e(1, 3)*cg1(3)*cg3(1) &
2321 +e(2, 1)*cg2(3)*cg1(1) &
2322 +e(2, 2)*cg2(3)*cg2(1) &
2323 +e(2, 3)*cg2(3)*cg3(1) &
2324 +e(3, 1)*cg3(3)*cg1(1) &
2325 +e(3, 2)*cg3(3)*cg2(1) )
2326
2327 !--------------------------------------------------
2328
2329 end do
2330 !--------------------------------------------------------------------
2331
2332 return
2333
2334 !####################################################################
2335 end subroutine elementstress_shell_mitc
2336 !####################################################################
2337
2338
2339 !####################################################################
2340 subroutine dl_shell &
2341 (etype, nn, ndof, xx, yy, zz, rho, thick, &
2342 ltype, params, vect, nsize, gausses)
2343 !####################################################################
2344
2345 use hecmw
2346 use m_utilities
2347 use mmechgauss
2349
2350 type(tgaussstatus), intent(in) :: gausses(:)
2351 !--------------------------------------------------------------------
2352
2353 integer(kind = kint), intent(in) :: etype
2354 integer(kind = kint), intent(in) :: nn
2355 integer(kind = kint), intent(in) :: ndof
2356 real(kind = kreal), intent(in) :: xx(*), yy(*), zz(*)
2357 real(kind = kreal), intent(in) :: rho
2358 real(kind = kreal), intent(in) :: thick
2359 real(kind = kreal), intent(in) :: params(*)
2360 real(kind = kreal), intent(out) :: vect(*)
2361 integer(kind = kint), intent(out) :: nsize
2362
2363 !--------------------------------------------------------------------
2364
2365 integer :: ivol, isurf
2366 integer :: lx, ly
2367 integer :: fetype
2368 integer :: ny
2369 integer :: i, m
2370 integer :: na, nb
2371 integer :: isize
2372 integer :: jsize1, jsize2, jsize3, &
2373 jsize4, jsize5, jsize6
2374 integer :: ltype
2375 integer :: n_totlyr, n_layer
2376
2377 real(kind = kreal) :: elem(3, nn)
2378 real(kind = kreal) :: val
2379 real(kind = kreal) :: ax, ay, az
2380 real(kind = kreal) :: rx, ry, rz
2381 real(kind = kreal) :: xi_lx, eta_lx, zeta_ly
2382 real(kind = kreal) :: w_w_lx, w_ly
2383 real(kind = kreal) :: naturalcoord(2)
2384 real(kind = kreal) :: nncoord(nn, 2)
2385 real(kind = kreal) :: shapefunc(nn)
2386 real(kind = kreal) :: shapederiv(nn, 2)
2387 real(kind = kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
2388 real(kind = kreal) :: v1_abs, v2_abs, v3_abs
2389 real(kind = kreal) :: a_over_2_v3(3, nn)
2390 real(kind = kreal) :: u_rot(3, nn)
2391 real(kind = kreal) :: dudxi_rot(3, nn), dudeta_rot(3, nn), &
2392 dudzeta_rot(3, nn)
2393 real(kind = kreal) :: g1(3), g2(3), g3(3)
2394 real(kind = kreal) :: g1_cross_g2(3)
2395 real(kind = kreal) :: e_0(3)
2396 real(kind = kreal) :: det
2397 real(kind = kreal) :: det_cg3(3)
2398 real(kind = kreal) :: det_cg3_abs
2399 real(kind = kreal) :: w_w_w_det
2400 real(kind = kreal) :: n(3, ndof*nn)
2401 real(kind = kreal) :: hx, hy, hz
2402 real(kind = kreal) :: phx, phy, phz
2403 real(kind = kreal) :: coefx, coefy, coefz
2404 real(kind = kreal) :: x, y, z
2405 real(kind = kreal) :: sumlyr
2406
2407 ny = 0
2408
2409 !--------------------------------------------------------------------
2410
2411 ! BX LTYPE=1 :BODY FORCE IN X-DIRECTION
2412 ! BY LTYPE=2 :BODY FORCE IN Y-DIRECTION
2413 ! BZ LTYPE=3 :BODY FORCE IN Z-DIRECTION
2414 ! CRAV LTYPE=4 :GRAVITY FORCE
2415 ! CENT LTYPE=5 :CENTRIFUGAL LOAD
2416 ! P LTYPE=10 :TRACTION IN NORMAL-DIRECTION FOR SHELL SURFACE
2417
2418 !--------------------------------------------------------------------
2419
2420 nsize = ndof*nn
2421
2422 !--------------------------------------------------------------------
2423
2424 val = params(1)
2425 ax = params(2)
2426 ay = params(3)
2427 az = params(4)
2428 rx = params(5)
2429 ry = params(6)
2430 rz = params(7)
2431
2432 !--------------------------------------------------------------------
2433
2434 ! MITC4
2435 if( etype .EQ. fe_mitc4_shell ) then
2436
2437 fetype = fe_mitc4_shell
2438
2439 ny = 2
2440
2441 ! MITC9
2442 else if( etype .EQ. fe_mitc9_shell ) then
2443
2444 fetype = fe_mitc9_shell
2445
2446 ny = 3
2447
2448 ! MITC3
2449 else if( etype .EQ. fe_mitc3_shell ) then
2450
2451 fetype = fe_mitc3_shell
2452
2453 ny = 2
2454
2455 end if
2456
2457 !--------------------------------------------------------------------
2458
2459 do na = 1, nn
2460
2461 elem(1, na) = xx(na)
2462 elem(2, na) = yy(na)
2463 elem(3, na) = zz(na)
2464
2465 end do
2466
2467 !-------------------------------------------------------------------
2468
2469 ! xi-coordinate at a node in a local element
2470 ! eta-coordinate at a node in a local element
2471 call getnodalnaturalcoord(fetype, nncoord)
2472
2473 !--------------------------------------------------------------------
2474
2475 ! Local load vector
2476 do isize = 1, ndof*nn
2477
2478 vect(isize) = 0.0d0
2479
2480 end do
2481
2482 !--------------------------------------------------------------------
2483
2484 ! xi-coordinate at the center point in a local element
2485 ! eta-coordinate at the center point in a local element
2486 naturalcoord(1) = 0.0d0
2487 naturalcoord(2) = 0.0d0
2488
2489 call getshapederiv(fetype, naturalcoord, shapederiv)
2490
2491 !--------------------------------------------------------------
2492
2493 ! Covariant basis vector
2494 do i = 1, 3
2495
2496 g1(i) = 0.0d0
2497
2498 do na = 1, nn
2499
2500 g1(i) = g1(i)+shapederiv(na, 1) &
2501 *elem(i, na)
2502
2503 end do
2504
2505 end do
2506
2507 e_0(1) = g1(1)
2508 e_0(2) = g1(2)
2509 e_0(3) = g1(3)
2510
2511 !--------------------------------------------------------------
2512
2513 do nb = 1, nn
2514
2515 !--------------------------------------------------------
2516
2517 naturalcoord(1) = nncoord(nb, 1)
2518 naturalcoord(2) = nncoord(nb, 2)
2519
2520 call getshapederiv(fetype, naturalcoord, shapederiv)
2521
2522 !--------------------------------------------------------
2523
2524 ! Covariant basis vector
2525 do i = 1, 3
2526
2527 g1(i) = 0.0d0
2528 g2(i) = 0.0d0
2529
2530 do na = 1, nn
2531
2532 g1(i) = g1(i)+shapederiv(na, 1) &
2533 *elem(i, na)
2534 g2(i) = g2(i)+shapederiv(na, 2) &
2535 *elem(i, na)
2536
2537 end do
2538
2539 end do
2540
2541 !--------------------------------------------------------
2542
2543 det_cg3(1) = g1(2)*g2(3)-g1(3)*g2(2)
2544 det_cg3(2) = g1(3)*g2(1)-g1(1)*g2(3)
2545 det_cg3(3) = g1(1)*g2(2)-g1(2)*g2(1)
2546
2547 det_cg3_abs = dsqrt( det_cg3(1)*det_cg3(1) &
2548 +det_cg3(2)*det_cg3(2) &
2549 +det_cg3(3)*det_cg3(3) )
2550
2551 v3(1, nb) = det_cg3(1)/det_cg3_abs
2552 v3(2, nb) = det_cg3(2)/det_cg3_abs
2553 v3(3, nb) = det_cg3(3)/det_cg3_abs
2554
2555 !--------------------------------------------------------
2556
2557 v2(1, nb) = v3(2, nb)*e_0(3)-v3(3, nb)*e_0(2)
2558 v2(2, nb) = v3(3, nb)*e_0(1)-v3(1, nb)*e_0(3)
2559 v2(3, nb) = v3(1, nb)*e_0(2)-v3(2, nb)*e_0(1)
2560
2561 v2_abs = dsqrt( v2(1, nb)*v2(1, nb) &
2562 +v2(2, nb)*v2(2, nb) &
2563 +v2(3, nb)*v2(3, nb) )
2564
2565 if( v2_abs .GT. 1.0d-15 ) then
2566
2567 v2(1, nb) = v2(1, nb)/v2_abs
2568 v2(2, nb) = v2(2, nb)/v2_abs
2569 v2(3, nb) = v2(3, nb)/v2_abs
2570
2571 v1(1, nb) = v2(2, nb)*v3(3, nb) &
2572 -v2(3, nb)*v3(2, nb)
2573 v1(2, nb) = v2(3, nb)*v3(1, nb) &
2574 -v2(1, nb)*v3(3, nb)
2575 v1(3, nb) = v2(1, nb)*v3(2, nb) &
2576 -v2(2, nb)*v3(1, nb)
2577
2578 v1_abs = dsqrt( v1(1, nb)*v1(1, nb) &
2579 +v1(2, nb)*v1(2, nb) &
2580 +v1(3, nb)*v1(3, nb) )
2581
2582 v1(1, nb) = v1(1, nb)/v1_abs
2583 v1(2, nb) = v1(2, nb)/v1_abs
2584 v1(3, nb) = v1(3, nb)/v1_abs
2585
2586 else
2587
2588 v1(1, nb) = 0.0d0
2589 v1(2, nb) = 0.0d0
2590 v1(3, nb) = -1.0d0
2591
2592 v2(1, nb) = 0.0d0
2593 v2(2, nb) = 1.0d0
2594 v2(3, nb) = 0.0d0
2595
2596 end if
2597
2598 !--------------------------------------------------------
2599
2600 v3(1, nb) = v1(2, nb)*v2(3, nb) &
2601 -v1(3, nb)*v2(2, nb)
2602 v3(2, nb) = v1(3, nb)*v2(1, nb) &
2603 -v1(1, nb)*v2(3, nb)
2604 v3(3, nb) = v1(1, nb)*v2(2, nb) &
2605 -v1(2, nb)*v2(1, nb)
2606
2607 v3_abs = dsqrt( v3(1, nb)*v3(1, nb) &
2608 +v3(2, nb)*v3(2, nb) &
2609 +v3(3, nb)*v3(3, nb) )
2610
2611 v3(1, nb) = v3(1, nb)/v3_abs
2612 v3(2, nb) = v3(2, nb)/v3_abs
2613 v3(3, nb) = v3(3, nb)/v3_abs
2614
2615 !--------------------------------------------------------
2616
2617 a_over_2_v3(1, nb) = 0.5d0*thick*v3(1, nb)
2618 a_over_2_v3(2, nb) = 0.5d0*thick*v3(2, nb)
2619 a_over_2_v3(3, nb) = 0.5d0*thick*v3(3, nb)
2620
2621 !--------------------------------------------------------
2622
2623 end do
2624
2625 !--------------------------------------------------------------------
2626
2627 ! Selction of load type
2628
2629 ivol = 0
2630 isurf = 0
2631
2632 if( ltype .LT. 10 ) then
2633
2634 ivol = 1
2635
2636 else if( ltype .GE. 10 ) then
2637
2638 isurf = 1
2639
2640 end if
2641
2642 !--------------------------------------------------------------------
2643
2644 !** Surface load
2645 if( isurf .EQ. 1 ) then
2646
2647 !--------------------------------------------------------
2648
2649 do lx = 1, numofquadpoints(fetype)
2650
2651 !--------------------------------------------------
2652
2653 call getquadpoint(fetype, lx, naturalcoord)
2654
2655 xi_lx = naturalcoord(1)
2656 eta_lx = naturalcoord(2)
2657
2658 w_w_lx = getweight(fetype, lx)
2659
2660 call getshapefunc(fetype, naturalcoord, shapefunc)
2661
2662 call getshapederiv(fetype, naturalcoord, shapederiv)
2663
2664 !--------------------------------------------------
2665
2666 do na = 1, nn
2667
2668 do i = 1, 3
2669
2670 u_rot(i, na) &
2671 = shapefunc(na) &
2672 *( 0.0d0*a_over_2_v3(i, na) )
2673
2674 dudxi_rot(i, na) &
2675 = shapederiv(na, 1) &
2676 *( 0.0d0*a_over_2_v3(i, na) )
2677 dudeta_rot(i, na) &
2678 = shapederiv(na, 2) &
2679 *( 0.0d0*a_over_2_v3(i, na) )
2680 dudzeta_rot(i, na) &
2681 = shapefunc(na) &
2682 *( a_over_2_v3(i, na) )
2683
2684 end do
2685
2686 end do
2687
2688 !--------------------------------------------------
2689
2690 ! Covariant basis vector
2691 do i = 1, 3
2692
2693 g1(i) = 0.0d0
2694 g2(i) = 0.0d0
2695 !g3(i) = 0.0D0
2696
2697 do na = 1, nn
2698
2699 g1(i) = g1(i)+shapederiv(na, 1) &
2700 *elem(i, na) &
2701 +dudxi_rot(i, na)
2702 g2(i) = g2(i)+shapederiv(na, 2) &
2703 *elem(i, na) &
2704 +dudeta_rot(i, na)
2705 !g3(i) = g3(i)+dudzeta_rot(i, na)
2706
2707 end do
2708
2709 end do
2710
2711 !--------------------------------------------------
2712
2713 !g3_abs = DSQRT( g3(1)*g3(1) &
2714 ! +g3(2)*g3(2) &
2715 ! +g3(3)*g3(3) )
2716
2717 !--------------------------------------------------
2718
2719 !e3_hat(1) = g3(1)/g3_abs
2720 !e3_hat(2) = g3(2)/g3_abs
2721 !e3_hat(3) = g3(3)/g3_abs
2722
2723 !--------------------------------------------------
2724
2725 ! Jacobian
2726 !det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
2727 ! +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
2728 ! +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )
2729
2730 !--------------------------------------------------
2731
2732 g1_cross_g2(1) = g1(2)*g2(3)-g1(3)*g2(2)
2733 g1_cross_g2(2) = g1(3)*g2(1)-g1(1)*g2(3)
2734 g1_cross_g2(3) = g1(1)*g2(2)-g1(2)*g2(1)
2735
2736 !--------------------------------------------------
2737
2738 do nb = 1, nn
2739
2740 jsize1 = ndof*(nb-1)+1
2741 jsize2 = ndof*(nb-1)+2
2742 jsize3 = ndof*(nb-1)+3
2743 jsize4 = ndof*(nb-1)+4
2744 jsize5 = ndof*(nb-1)+5
2745 jsize6 = ndof*(nb-1)+6
2746
2747 n(1, jsize1) = shapefunc(nb)
2748 n(1, jsize2) = 0.0d0
2749 n(1, jsize3) = 0.0d0
2750 n(1, jsize4) = 0.0d0
2751 n(1, jsize5) = 0.0d0
2752 n(1, jsize6) = 0.0d0
2753 n(2, jsize1) = 0.0d0
2754 n(2, jsize2) = shapefunc(nb)
2755 n(2, jsize3) = 0.0d0
2756 n(2, jsize4) = 0.0d0
2757 n(2, jsize5) = 0.0d0
2758 n(2, jsize6) = 0.0d0
2759 n(3, jsize1) = 0.0d0
2760 n(3, jsize2) = 0.0d0
2761 n(3, jsize3) = shapefunc(nb)
2762 n(3, jsize4) = 0.0d0
2763 n(3, jsize5) = 0.0d0
2764 n(3, jsize6) = 0.0d0
2765
2766 end do
2767
2768 do isize = 1, ndof*nn
2769
2770 vect(isize) &
2771 = vect(isize) &
2772 +w_w_lx*( n(1, isize)*g1_cross_g2(1) &
2773 +n(2, isize)*g1_cross_g2(2) &
2774 +n(3, isize)*g1_cross_g2(3) )*val
2775
2776 end do
2777
2778 !--------------------------------------------------
2779
2780 end do
2781
2782 !--------------------------------------------------------
2783
2784 end if
2785
2786 !--------------------------------------------------------------------
2787
2788 !** Volume load
2789 if( ivol .EQ. 1 ) then
2790
2791 !--------------------------------------------------------
2792 n_totlyr = gausses(1)%pMaterial%totallyr
2793 do n_layer=1,n_totlyr
2794 do ly = 1, ny
2795
2796 !--------------------------------------------------
2797
2798
2799 sumlyr = 0.0d0
2800 do m = 1, n_layer
2801 sumlyr = sumlyr + 2 * gausses(1)%pMaterial%shell_var(m)%weight
2802 end do
2803 zeta_ly = -1 + sumlyr - gausses(1)%pMaterial%shell_var(n_layer)%weight*(1-xg(ny, ly))
2804
2805 !zeta_ly = xg(ny, ly)
2806 w_ly = wgt(ny, ly)
2807
2808 !--------------------------------------------------
2809
2810 do lx = 1, numofquadpoints(fetype)
2811
2812 !--------------------------------------------
2813
2814 call getquadpoint(fetype, lx, naturalcoord)
2815
2816 xi_lx = naturalcoord(1)
2817 eta_lx = naturalcoord(2)
2818
2819 w_w_lx = getweight(fetype, lx)
2820
2821 call getshapefunc(fetype, naturalcoord, shapefunc)
2822
2823 call getshapederiv(fetype, naturalcoord, shapederiv)
2824
2825 !--------------------------------------------
2826
2827 do na = 1, nn
2828
2829 do i = 1, 3
2830
2831 u_rot(i, na) &
2832 = shapefunc(na) &
2833 *( zeta_ly*a_over_2_v3(i, na) )
2834
2835 dudxi_rot(i, na) &
2836 = shapederiv(na, 1) &
2837 *( zeta_ly*a_over_2_v3(i, na) )
2838 dudeta_rot(i, na) &
2839 = shapederiv(na, 2) &
2840 *( zeta_ly*a_over_2_v3(i, na) )
2841 dudzeta_rot(i, na) &
2842 = shapefunc(na) &
2843 *( a_over_2_v3(i, na) )
2844
2845 end do
2846
2847 end do
2848
2849 !--------------------------------------------
2850
2851 ! Covariant basis vector
2852 do i = 1, 3
2853
2854 g1(i) = 0.0d0
2855 g2(i) = 0.0d0
2856 g3(i) = 0.0d0
2857
2858 do na = 1, nn
2859
2860 g1(i) = g1(i)+shapederiv(na, 1) &
2861 *elem(i, na) &
2862 +dudxi_rot(i, na)
2863 g2(i) = g2(i)+shapederiv(na, 2) &
2864 *elem(i, na) &
2865 +dudeta_rot(i, na)
2866 g3(i) = g3(i)+dudzeta_rot(i, na)
2867
2868 end do
2869
2870 end do
2871
2872 !--------------------------------------------
2873
2874 ! Jacobian
2875 det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
2876 +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
2877 +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )
2878
2879 !--------------------------------------------
2880
2881 ! [ N ] matrix
2882 do nb = 1, nn
2883
2884 jsize1 = ndof*(nb-1)+1
2885 jsize2 = ndof*(nb-1)+2
2886 jsize3 = ndof*(nb-1)+3
2887 jsize4 = ndof*(nb-1)+4
2888 jsize5 = ndof*(nb-1)+5
2889 jsize6 = ndof*(nb-1)+6
2890
2891 n(1, jsize1) = shapefunc(nb)
2892 n(2, jsize1) = 0.0d0
2893 n(3, jsize1) = 0.0d0
2894 n(1, jsize2) = 0.0d0
2895 n(2, jsize2) = shapefunc(nb)
2896 n(3, jsize2) = 0.0d0
2897 n(1, jsize3) = 0.0d0
2898 n(2, jsize3) = 0.0d0
2899 n(3, jsize3) = shapefunc(nb)
2900 n(1, jsize4) = 0.0d0
2901 n(2, jsize4) = -u_rot(3, nb)
2902 n(3, jsize4) = u_rot(2, nb)
2903 n(1, jsize5) = u_rot(3, nb)
2904 n(2, jsize5) = 0.0d0
2905 n(3, jsize5) = -u_rot(1, nb)
2906 n(1, jsize6) = -u_rot(2, nb)
2907 n(2, jsize6) = u_rot(1, nb)
2908 n(3, jsize6) = 0.0d0
2909
2910 enddo
2911
2912 !--------------------------------------------
2913
2914 w_w_w_det = w_w_lx*w_ly*det
2915
2916 !--------------------------------------------
2917
2918 if( ltype .EQ. 1 ) then
2919
2920 do isize = 1, ndof*nn
2921
2922 vect(isize) = vect(isize)+w_w_w_det*n(1, isize)*val
2923
2924 end do
2925
2926 else if( ltype .EQ. 2 ) then
2927
2928 do isize = 1, ndof*nn
2929
2930 vect(isize) = vect(isize)+w_w_w_det*n(2, isize)*val
2931
2932 end do
2933
2934 else if( ltype .EQ. 3 ) then
2935
2936 do isize = 1, ndof*nn
2937
2938 vect(isize) = vect(isize)+w_w_w_det*n(3, isize)*val
2939
2940 end do
2941
2942 else if( ltype .EQ. 4 ) then
2943
2944 do isize = 1, ndof*nn
2945
2946 vect(isize) = vect(isize)+w_w_w_det*rho*ax*n(1, isize)*val
2947 vect(isize) = vect(isize)+w_w_w_det*rho*ay*n(2, isize)*val
2948 vect(isize) = vect(isize)+w_w_w_det*rho*az*n(3, isize)*val
2949
2950 end do
2951
2952 else if( ltype .EQ. 5 ) then
2953
2954 x = 0.0d0
2955 y = 0.0d0
2956 z = 0.0d0
2957
2958 do nb = 1, nn
2959
2960 x = x+shapefunc(nb)*elem(1, nb)
2961 y = y+shapefunc(nb)*elem(2, nb)
2962 z = z+shapefunc(nb)*elem(3, nb)
2963
2964 end do
2965
2966 hx = ax+( (x-ax)*rx+(y-ay)*ry+(z-az)*rz )/( rx**2+ry**2+rz**2 )*rx
2967 hy = ay+( (x-ax)*rx+(y-ay)*ry+(z-az)*rz )/( rx**2+ry**2+rz**2 )*ry
2968 hz = az+( (x-ax)*rx+(y-ay)*ry+(z-az)*rz )/( rx**2+ry**2+rz**2 )*rz
2969
2970 phx = x-hx
2971 phy = y-hy
2972 phz = z-hz
2973
2974 coefx = phx*val*rho*val
2975 coefy = phy*val*rho*val
2976 coefz = phz*val*rho*val
2977
2978 do isize = 1, ndof*nn
2979
2980 vect(isize) &
2981 = vect(isize) &
2982 +w_w_w_det*( n(1, isize)*coefx &
2983 +n(2, isize)*coefy &
2984 +n(3, isize)*coefz )
2985
2986 end do
2987
2988 end if
2989
2990 !--------------------------------------------
2991
2992 end do
2993
2994 !----------------------------------------------
2995
2996 end do
2997
2998 !----------------------------------------------
2999
3000 end do
3001
3002 !--------------------------------------------------------
3003
3004 end if
3005 !--------------------------------------------------------------------
3006
3007 return
3008
3009 !####################################################################
3010 end subroutine dl_shell
3011 !####################################################################
3012
3013
3014 !####################################################################
3015 subroutine dl_shell_33 &
3016 (ic_type, nn, ndof, xx, yy, zz, rho, thick, &
3017 ltype, params, vect, nsize, gausses)
3018 !####################################################################
3019
3020 use hecmw
3021 use m_utilities
3022 use mmechgauss
3024
3025 type(tgaussstatus) :: gausses(:)
3026 !--------------------------------------------------------------------
3027
3028 integer(kind = kint) :: ic_type
3029 integer(kind = kint) :: nn
3030 integer(kind = kint) :: ndof
3031 real(kind = kreal) :: xx(*), yy(*), zz(*)
3032 real(kind = kreal) :: rho
3033 real(kind = kreal) :: thick
3034 real(kind = kreal) :: params(*)
3035 real(kind = kreal) :: vect(*)
3036 integer(kind = kint) :: nsize
3037 integer :: ltype, i
3038 real(kind = kreal) :: tmp(24)
3039 !--------------------------------------------------------------------
3040
3041 if(ic_type == 761)then
3042 !ic_type = 731
3043 !nn = 3
3044 !ndof = 6
3045 call dl_shell(731, 3, 6, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
3046 !ic_type = 761
3047 !nn = 6
3048 !ndof = 3
3049
3050 tmp = 0.0
3051 do i=1,18
3052 tmp(i) = vect(i)
3053 enddo
3054
3055 vect( 1) = tmp(1)
3056 vect( 2) = tmp(2)
3057 vect( 3) = tmp(3)
3058 vect( 4) = tmp(7)
3059 vect( 5) = tmp(8)
3060 vect( 6) = tmp(9)
3061 vect( 7) = tmp(13)
3062 vect( 8) = tmp(14)
3063 vect( 9) = tmp(15)
3064 vect(10) = tmp(4)
3065 vect(11) = tmp(5)
3066 vect(12) = tmp(6)
3067 vect(13) = tmp(10)
3068 vect(14) = tmp(11)
3069 vect(15) = tmp(12)
3070 vect(16) = tmp(16)
3071 vect(17) = tmp(17)
3072 vect(18) = tmp(18)
3073
3074 elseif(ic_type == 781)then
3075 !ic_type = 741
3076 !nn = 4
3077 !ndof = 6
3078 call dl_shell(741, 4, 6, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
3079 !ic_type = 781
3080 !nn = 8
3081 !ndof = 3
3082
3083 tmp = 0.0
3084 do i=1,24
3085 tmp(i) = vect(i)
3086 enddo
3087
3088 vect( 1) = tmp(1)
3089 vect( 2) = tmp(2)
3090 vect( 3) = tmp(3)
3091 vect( 4) = tmp(7)
3092 vect( 5) = tmp(8)
3093 vect( 6) = tmp(9)
3094 vect( 7) = tmp(13)
3095 vect( 8) = tmp(14)
3096 vect( 9) = tmp(15)
3097 vect(10) = tmp(19)
3098 vect(11) = tmp(20)
3099 vect(12) = tmp(21)
3100 vect(13) = tmp(4)
3101 vect(14) = tmp(5)
3102 vect(15) = tmp(6)
3103 vect(16) = tmp(10)
3104 vect(17) = tmp(11)
3105 vect(18) = tmp(12)
3106 vect(19) = tmp(16)
3107 vect(20) = tmp(17)
3108 vect(21) = tmp(18)
3109 vect(22) = tmp(22)
3110 vect(23) = tmp(23)
3111 vect(24) = tmp(24)
3112
3113 endif
3114
3115 end subroutine dl_shell_33
3116
3117 !####################################################################
3118 ! this subroutine can be used only for linear analysis
3120 (etype, nn, ndof, ecoord, u, du, gausses, qf, thick, mixflag, nddisp)
3121 !####################################################################
3122
3123 use mmechgauss
3125 use m_matmatrix
3126
3127 !--------------------------------------------------------------------
3128
3129 integer(kind = kint), intent(in) :: etype
3130 integer(kind = kint), intent(in) :: nn, mixflag
3131 integer(kind = kint), intent(in) :: ndof
3132 real(kind = kreal), intent(in) :: ecoord(3, nn)
3133 real(kind = kreal), intent(in) :: u(:, :)
3134 real(kind = kreal), intent(in) :: du(:, :)
3135 type(tgaussstatus), intent(in) :: gausses(:)
3136 real(kind = kreal), intent(out) :: qf(:)
3137 real(kind = kreal), intent(in) :: thick
3138
3139 real(kind = kreal), intent(in), optional :: nddisp(3, nn)
3140 !--------------------------------------------------------------------
3141
3142 real(kind = kreal) :: stiff(nn*ndof, nn*ndof), totaldisp(nn*ndof)
3143 integer(kind = kint) :: i
3144
3145 call stf_shell_mitc(etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag, nddisp)
3146
3147 totaldisp = 0.d0
3148 do i=1,nn
3149 totaldisp(ndof*(i-1)+1:ndof*i) = u(1:ndof,i) + du(1:ndof,i)
3150 end do
3151
3152 qf = matmul(stiff,totaldisp)
3153
3154 end subroutine updatest_shell_mitc
3155
3156 !####################################################################
3157 ! this subroutine can be used only for linear analysis
3159 (etype, nn, ndof, ecoord, u, du, gausses, qf, thick, mixflag, nddisp)
3160 !####################################################################
3161
3162 use mmechgauss
3164 use m_matmatrix
3165
3166 !--------------------------------------------------------------------
3167
3168 integer(kind = kint), intent(in) :: etype
3169 integer(kind = kint), intent(in) :: nn, mixflag
3170 integer(kind = kint), intent(in) :: ndof
3171 real(kind = kreal), intent(in) :: ecoord(3, nn)
3172 real(kind = kreal), intent(in) :: u(3, nn*2)
3173 real(kind = kreal), intent(in) :: du(3, nn*2)
3174 type(tgaussstatus), intent(in) :: gausses(:)
3175 real(kind = kreal), intent(out) :: qf(:)
3176 real(kind = kreal), intent(in) :: thick
3177
3178 real(kind = kreal), intent(in), optional :: nddisp(3, nn)
3179 !--------------------------------------------------------------------
3180
3181 real(kind = kreal) :: stiff(nn*ndof, nn*ndof), totaldisp(nn*ndof)
3182 integer(kind = kint) :: i
3183
3184 call stf_shell_mitc(etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag, nddisp)
3185
3186 totaldisp = 0.d0
3187 do i=1,nn
3188 totaldisp(ndof*(i-1)+1:ndof*(i-1)+3) = u(1:3,2*i-1) + du(1:3,2*i-1)
3189 totaldisp(ndof*(i-1)+4:ndof*(i-1)+6) = u(1:3,2*i) + du(1:3,2*i)
3190 end do
3191
3192 qf = matmul(stiff,totaldisp)
3193
3194 end subroutine updatest_shell_mitc33
3195
3196 !####################################################################
3197 subroutine mass_shell(etype, nn, elem, rho, thick, gausses, mass, lumped)
3198 !####################################################################
3199 use hecmw
3200 use m_utilities
3201 use mmechgauss
3203 integer(kind = kint), intent(in) :: etype
3204 integer(kind = kint), intent(in) :: nn
3205 real(kind = kreal), intent(in) :: elem(3,nn)
3206 real(kind = kreal), intent(in) :: rho
3207 real(kind = kreal), intent(in) :: thick
3208 type(tgaussstatus), intent(in) :: gausses(:)
3209 real(kind=kreal), intent(out) :: mass(:,:)
3210 real(kind=kreal), intent(out) :: lumped(:)
3211
3212 !--------------------------------------------------------------------
3213
3214 integer :: lx, ly, nsize, ndof
3215 integer :: fetype
3216 integer :: ny
3217 integer :: i, m
3218 integer :: na, nb
3219 integer :: jsize1, jsize2, jsize3, jsize4, jsize5, jsize6
3220 integer :: n_totlyr, n_layer
3221
3222 real(kind = kreal) :: xi_lx, eta_lx, zeta_ly
3223 real(kind = kreal) :: w_w_lx, w_ly
3224 real(kind = kreal) :: naturalcoord(2)
3225 real(kind = kreal) :: nncoord(nn, 2)
3226 real(kind = kreal) :: shapefunc(nn)
3227 real(kind = kreal) :: shapederiv(nn, 2)
3228 real(kind = kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
3229 real(kind = kreal) :: v1_abs, v2_abs, v3_abs
3230 real(kind = kreal) :: a_over_2_v3(3, nn)
3231 real(kind = kreal) :: u_rot(3, nn)
3232 real(kind = kreal) :: dudxi_rot(3, nn), dudeta_rot(3, nn), dudzeta_rot(3, nn)
3233 real(kind = kreal) :: g1(3), g2(3), g3(3)
3234 real(kind = kreal) :: e_0(3)
3235 real(kind = kreal) :: det
3236 real(kind = kreal) :: det_cg3(3)
3237 real(kind = kreal) :: det_cg3_abs
3238 real(kind = kreal) :: w_w_w_det
3239 real(kind = kreal) :: n(3, 6*nn)
3240 real(kind = kreal) :: sumlyr, totalmass, totdiag
3241
3242 ny = 0; ndof=6
3243 nsize = ndof*nn
3244
3245 !--------------------------------------------------------------------
3246
3247 ! MITC4
3248 if( etype == fe_mitc4_shell ) then
3249 fetype = fe_mitc4_shell
3250 ny = 2
3251
3252 ! MITC9
3253 else if( etype == fe_mitc9_shell ) then
3254 fetype = fe_mitc9_shell
3255 ny = 3
3256
3257 ! MITC3
3258 else if( etype == fe_mitc3_shell ) then
3259 fetype = fe_mitc3_shell
3260 ny = 2
3261
3262 end if
3263
3264 !-------------------------------------------------------------------
3265
3266 ! xi-coordinate at a node in a local element
3267 ! eta-coordinate at a node in a local element
3268 call getnodalnaturalcoord(fetype, nncoord)
3269
3270 ! xi-coordinate at the center point in a local element
3271 ! eta-coordinate at the center point in a local element
3272 naturalcoord(1) = 0.0d0
3273 naturalcoord(2) = 0.0d0
3274
3275 call getshapederiv(fetype, naturalcoord, shapederiv)
3276
3277 !--------------------------------------------------------------
3278
3279 ! Covariant basis vector
3280 g1(:) = matmul( elem, shapederiv(:,1) )
3281
3282 e_0(:) = g1(:)
3283
3284 !--------------------------------------------------------------
3285
3286 do nb = 1, nn
3287
3288 !--------------------------------------------------------
3289
3290 naturalcoord(1) = nncoord(nb, 1)
3291 naturalcoord(2) = nncoord(nb, 2)
3292 call getshapederiv(fetype, naturalcoord, shapederiv)
3293
3294 !--------------------------------------------------------
3295
3296 ! Covariant basis vector
3297 g1(:) = matmul( elem, shapederiv(:,1) )
3298 g2(:) = matmul( elem, shapederiv(:,2) )
3299
3300 !--------------------------------------------------------
3301
3302 det_cg3(1) = g1(2)*g2(3)-g1(3)*g2(2)
3303 det_cg3(2) = g1(3)*g2(1)-g1(1)*g2(3)
3304 det_cg3(3) = g1(1)*g2(2)-g1(2)*g2(1)
3305
3306 det_cg3_abs = dsqrt( det_cg3(1)*det_cg3(1) &
3307 +det_cg3(2)*det_cg3(2) &
3308 +det_cg3(3)*det_cg3(3) )
3309
3310 v3(:, nb) = det_cg3(:)/det_cg3_abs
3311
3312 !--------------------------------------------------------
3313
3314 v2(1, nb) = v3(2, nb)*e_0(3)-v3(3, nb)*e_0(2)
3315 v2(2, nb) = v3(3, nb)*e_0(1)-v3(1, nb)*e_0(3)
3316 v2(3, nb) = v3(1, nb)*e_0(2)-v3(2, nb)*e_0(1)
3317
3318 v2_abs = dsqrt( v2(1, nb)*v2(1, nb) &
3319 +v2(2, nb)*v2(2, nb) &
3320 +v2(3, nb)*v2(3, nb) )
3321
3322 if( v2_abs > 1.0d-15 ) then
3323
3324 v2(1, nb) = v2(1, nb)/v2_abs
3325 v2(2, nb) = v2(2, nb)/v2_abs
3326 v2(3, nb) = v2(3, nb)/v2_abs
3327
3328 v1(1, nb) = v2(2, nb)*v3(3, nb) &
3329 -v2(3, nb)*v3(2, nb)
3330 v1(2, nb) = v2(3, nb)*v3(1, nb) &
3331 -v2(1, nb)*v3(3, nb)
3332 v1(3, nb) = v2(1, nb)*v3(2, nb) &
3333 -v2(2, nb)*v3(1, nb)
3334
3335 v1_abs = dsqrt( v1(1, nb)*v1(1, nb) &
3336 +v1(2, nb)*v1(2, nb) &
3337 +v1(3, nb)*v1(3, nb) )
3338
3339 v1(1, nb) = v1(1, nb)/v1_abs
3340 v1(2, nb) = v1(2, nb)/v1_abs
3341 v1(3, nb) = v1(3, nb)/v1_abs
3342
3343 else
3344
3345 v1(1, nb) = 0.0d0
3346 v1(2, nb) = 0.0d0
3347 v1(3, nb) = -1.0d0
3348
3349 v2(1, nb) = 0.0d0
3350 v2(2, nb) = 1.0d0
3351 v2(3, nb) = 0.0d0
3352
3353 end if
3354
3355 !--------------------------------------------------------
3356
3357 v3(1, nb) = v1(2, nb)*v2(3, nb) &
3358 -v1(3, nb)*v2(2, nb)
3359 v3(2, nb) = v1(3, nb)*v2(1, nb) &
3360 -v1(1, nb)*v2(3, nb)
3361 v3(3, nb) = v1(1, nb)*v2(2, nb) &
3362 -v1(2, nb)*v2(1, nb)
3363
3364 v3_abs = dsqrt( v3(1, nb)*v3(1, nb) &
3365 +v3(2, nb)*v3(2, nb) &
3366 +v3(3, nb)*v3(3, nb) )
3367
3368 v3(1, nb) = v3(1, nb)/v3_abs
3369 v3(2, nb) = v3(2, nb)/v3_abs
3370 v3(3, nb) = v3(3, nb)/v3_abs
3371
3372 !--------------------------------------------------------
3373
3374 a_over_2_v3(1, nb) = 0.5d0*thick*v3(1, nb)
3375 a_over_2_v3(2, nb) = 0.5d0*thick*v3(2, nb)
3376 a_over_2_v3(3, nb) = 0.5d0*thick*v3(3, nb)
3377
3378 !--------------------------------------------------------
3379
3380 end do
3381
3382 !--------------------------------------------------------------------
3383
3384 mass(:,:) = 0.0d0
3385 totalmass = 0.d0
3386 n_totlyr = gausses(1)%pMaterial%totallyr
3387 do n_layer=1,n_totlyr
3388 do ly = 1, ny
3389
3390 !--------------------------------------------------
3391
3392 sumlyr = 0.0d0
3393 do m = 1, n_layer
3394 sumlyr = sumlyr + 2 * gausses(1)%pMaterial%shell_var(m)%weight
3395 end do
3396 zeta_ly = -1 + sumlyr - gausses(1)%pMaterial%shell_var(n_layer)%weight*(1-xg(ny, ly))
3397
3398 !zeta_ly = xg(ny, ly)
3399 w_ly = wgt(ny, ly)
3400
3401 !--------------------------------------------------
3402
3403 do lx = 1, numofquadpoints(fetype)
3404
3405 !--------------------------------------------
3406
3407 call getquadpoint(fetype, lx, naturalcoord)
3408
3409 xi_lx = naturalcoord(1)
3410 eta_lx = naturalcoord(2)
3411
3412 w_w_lx = getweight(fetype, lx)
3413
3414 call getshapefunc(fetype, naturalcoord, shapefunc)
3415 call getshapederiv(fetype, naturalcoord, shapederiv)
3416
3417 !--------------------------------------------
3418
3419 do na = 1, nn
3420
3421 do i = 1, 3
3422
3423 u_rot(i, na) = shapefunc(na)*( zeta_ly*a_over_2_v3(i, na) )
3424
3425 dudxi_rot(i, na) = shapederiv(na, 1) &
3426 *( zeta_ly*a_over_2_v3(i, na) )
3427 dudeta_rot(i, na) = shapederiv(na, 2) &
3428 *( zeta_ly*a_over_2_v3(i, na) )
3429 dudzeta_rot(i, na) = shapefunc(na) &
3430 *( a_over_2_v3(i, na) )
3431
3432 end do
3433
3434 end do
3435
3436 !--------------------------------------------
3437
3438 ! Covariant basis vector
3439 do i = 1, 3
3440 g1(i) = 0.0d0
3441 g2(i) = 0.0d0
3442 g3(i) = 0.0d0
3443 do na = 1, nn
3444 g1(i) = g1(i)+shapederiv(na, 1) *elem(i, na) &
3445 +dudxi_rot(i, na)
3446 g2(i) = g2(i)+shapederiv(na, 2) *elem(i, na) &
3447 +dudeta_rot(i, na)
3448 g3(i) = g3(i)+dudzeta_rot(i, na)
3449 end do
3450 end do
3451
3452 !--------------------------------------------
3453
3454 ! Jacobian
3455 det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
3456 +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
3457 +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )
3458
3459 !--------------------------------------------
3460
3461 ! [ N ] matrix
3462 do nb = 1, nn
3463
3464 jsize1 = ndof*(nb-1)+1
3465 jsize2 = ndof*(nb-1)+2
3466 jsize3 = ndof*(nb-1)+3
3467 jsize4 = ndof*(nb-1)+4
3468 jsize5 = ndof*(nb-1)+5
3469 jsize6 = ndof*(nb-1)+6
3470
3471 n(1, jsize1) = shapefunc(nb)
3472 n(2, jsize1) = 0.0d0
3473 n(3, jsize1) = 0.0d0
3474 n(1, jsize2) = 0.0d0
3475 n(2, jsize2) = shapefunc(nb)
3476 n(3, jsize2) = 0.0d0
3477 n(1, jsize3) = 0.0d0
3478 n(2, jsize3) = 0.0d0
3479 n(3, jsize3) = shapefunc(nb)
3480 n(1, jsize4) = 0.0d0
3481 n(2, jsize4) = -u_rot(3, nb)
3482 n(3, jsize4) = u_rot(2, nb)
3483 n(1, jsize5) = u_rot(3, nb)
3484 n(2, jsize5) = 0.0d0
3485 n(3, jsize5) = -u_rot(1, nb)
3486 n(1, jsize6) = -u_rot(2, nb)
3487 n(2, jsize6) = u_rot(1, nb)
3488 n(3, jsize6) = 0.0d0
3489
3490 enddo
3491
3492 !--------------------------------------------
3493
3494 w_w_w_det = w_w_lx*w_ly*det
3495 mass(1:nsize,1:nsize) = mass(1:nsize,1:nsize)+ matmul( transpose(n), n )*w_w_w_det*rho
3496 totalmass = totalmass + w_w_w_det*rho
3497 !--------------------------------------------
3498
3499 end do
3500
3501 !----------------------------------------------
3502
3503 end do
3504
3505 !----------------------------------------------
3506
3507 end do
3508 totalmass = totalmass*3.d0
3509
3510 totdiag=0.d0
3511 do nb = 1, nn
3512 DO i = 1, 6
3513 lx = (nb-1)*ndof+i
3514 totdiag = totdiag + mass(lx,lx)
3515 END DO
3516 ENDDO
3517 lumped(:)=0.d0
3518 do nb = 1, nn
3519 DO i = 1, 6
3520 lx = (nb-1)*ndof+i
3521 lumped(lx) = mass(lx,lx)/totdiag* totalmass
3522 END DO
3523 ENDDO
3524
3525 !####################################################################
3526 end subroutine mass_shell
3527 !####################################################################
3528
3529end module m_static_lib_shell
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
Definition: element.f90:640
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:482
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate deivatives of shape fucntion in natural coordiante system.
Definition: element.f90:571
integer, parameter fe_mitc4_shell
Definition: element.f90:92
integer, parameter fe_mitc9_shell
Definition: element.f90:94
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:528
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:445
integer, parameter fe_mitc3_shell
Definition: element.f90:91
subroutine getnodalnaturalcoord(fetype, nncoord)
Definition: element.f90:692
This module provides data for gauss quadrature.
Definition: GaussM.f90:6
real(kind=kreal), dimension(3, 3) wgt
wieght of gauss points
Definition: GaussM.f90:10
real(kind=kreal), dimension(3, 3) xg
abscissa of gauss points
Definition: GaussM.f90:9
Definition: hecmw.f90:6
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
subroutine matlmatrix_shell(gauss, secttype, d, e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, alpha, n_layer)
subroutine stf_shell_mitc(etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag, nddisp)
subroutine elementstress_shell_mitc(etype, nn, ndof, ecoord, gausses, edisp, strain, stress, thick, zeta, n_layer, n_totlyr)
subroutine dl_shell_33(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
subroutine updatest_shell_mitc33(etype, nn, ndof, ecoord, u, du, gausses, qf, thick, mixflag, nddisp)
subroutine mass_shell(etype, nn, elem, rho, thick, gausses, mass, lumped)
subroutine dl_shell(etype, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
subroutine updatest_shell_mitc(etype, nn, ndof, ecoord, u, du, gausses, qf, thick, mixflag, nddisp)
This module provides aux functions.
Definition: utilities.f90:6
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13