6 use hecmw,
only : kint, kreal
57 (etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag, nddisp)
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)
71 real(kind = kreal),
intent(out) :: stiff(:, :)
72 real(kind = kreal),
intent(in) :: thick
74 real(kind = kreal),
intent(in),
optional :: nddisp(3, nn)
78 integer :: flag, flag_dof
84 integer :: npoints_tying(3)
87 integer :: isize, jsize
88 integer :: jsize1, jsize2, jsize3, &
89 jsize4, jsize5, jsize6
90 integer :: n_layer,n_totlyr, sstable(24)
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), &
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), &
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
180 if(
present( nddisp ) )
then
182 unode(:, 1:nn) = nddisp(:, :)
188 flag = gausses(1)%pMaterial%nlgeom_flag
190 if( .not.
present( nddisp ) ) flag = infinitesimal
194 elem(:, :) = ecoord(:, :)
198 tmpstiff(:, :) = 0.0d0
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
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
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 )
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 )
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 )
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 )
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 )
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 )
275 xxi_di(1, 1) = -1.0d0
278 xxi_di(4, 1) = -1.0d0
280 xxi_di(6, 1) = -1.0d0
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
290 xxi_di(1, 2) = -1.0d0
295 xxi_di(6, 2) = -1.0d0
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
312 tpcoord(1, 1, 1) = 0.5d0
313 tpcoord(2, 1, 1) = 0.0d0
314 tpcoord(3, 1, 1) = 0.5d0
316 tpcoord(1, 2, 1) = 0.0d0
317 tpcoord(2, 2, 1) = 0.5d0
318 tpcoord(3, 2, 1) = 0.5d0
328 naturalcoord(1) = 0.0d0
329 naturalcoord(2) = 0.0d0
342 g1(i) = g1(i)+shapederiv(na, 1) &
359 naturalcoord(1) = nncoord(nb, 1)
360 naturalcoord(2) = nncoord(nb, 2)
374 g1(i) = g1(i)+shapederiv(na, 1) &
376 g2(i) = g2(i)+shapederiv(na, 2) &
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)
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) )
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
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)
403 v2_abs = dsqrt( v2(1, nb)*v2(1, nb) &
404 +v2(2, nb)*v2(2, nb) &
405 +v2(3, nb)*v2(3, nb) )
407 if( v2_abs .GT. 1.0d-15 )
then
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
413 v1(1, nb) = v2(2, nb)*v3(3, nb) &
415 v1(2, nb) = v2(3, nb)*v3(1, nb) &
417 v1(3, nb) = v2(1, nb)*v3(2, nb) &
420 v1_abs = dsqrt( v1(1, nb)*v1(1, nb) &
421 +v1(2, nb)*v1(2, nb) &
422 +v1(3, nb)*v1(3, nb) )
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
442 v3(1, nb) = v1(2, nb)*v2(3, nb) &
444 v3(2, nb) = v1(3, nb)*v2(1, nb) &
446 v3(3, nb) = v1(1, nb)*v2(2, nb) &
449 v3_abs = dsqrt( v3(1, nb)*v3(1, nb) &
450 +v3(2, nb)*v3(2, nb) &
451 +v3(3, nb)*v3(3, nb) )
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
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)
472 n_totlyr = gausses(1)%pMaterial%totallyr
473 do n_layer=1,n_totlyr
499 do ip = 1, npoints_tying(it)
503 naturalcoord(1) = tpcoord(ip, 1, it)
504 naturalcoord(2) = tpcoord(ip, 2, it)
518 *( zeta_ly*a_over_2_v3(i, na) )
521 = shapederiv(na, 1) &
522 *( zeta_ly*a_over_2_v3(i, na) )
524 = shapederiv(na, 2) &
525 *( zeta_ly*a_over_2_v3(i, na) )
528 *( a_over_2_v3(i, na) )
545 g1(i) = g1(i)+shapederiv(na, 1) &
548 g2(i) = g2(i)+shapederiv(na, 2) &
551 g3(i) = g3(i)+dudzeta_rot(i, na)
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
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
652 sumlyr = sumlyr +2 *gausses(1)%pMaterial%shell_var(m)%weight
654 zeta_ly = -1 +sumlyr -gausses(1)%pMaterial%shell_var(n_layer)%weight*(1-
xg(ny, ly))
666 xi_lx = naturalcoord(1)
667 eta_lx = naturalcoord(2)
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)
701 *( zeta_ly*a_over_2_v3(i, na) )
704 = shapederiv(na, 1) &
705 *( zeta_ly*a_over_2_v3(i, na) )
707 = shapederiv(na, 2) &
708 *( zeta_ly*a_over_2_v3(i, na) )
711 *( a_over_2_v3(i, na) )
728 g1(i) = g1(i)+shapederiv(na, 1) &
731 g2(i) = g2(i)+shapederiv(na, 2) &
734 g3(i) = g3(i)+dudzeta_rot(i, na)
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) )
753 *( g2(2)*g3(3)-g2(3)*g3(2) )
755 *( g2(3)*g3(1)-g2(1)*g3(3) )
757 *( g2(1)*g3(2)-g2(2)*g3(1) )
759 *( g3(2)*g1(3)-g3(3)*g1(2) )
761 *( g3(3)*g1(1)-g3(1)*g1(3) )
763 *( g3(1)*g1(2)-g3(2)*g1(1) )
765 *( g1(2)*g2(3)-g1(3)*g2(2) )
767 *( g1(3)*g2(1)-g1(1)*g2(3) )
769 *( g1(1)*g2(2)-g1(2)*g2(1) )
773 g3_abs = dsqrt( g3(1)*g3(1) &
781 e3_hat(1) = g3(1)/g3_abs
782 e3_hat(2) = g3(2)/g3_abs
783 e3_hat(3) = g3(3)/g3_abs
785 e1_hat(1) = g2(2)*e3_hat(3) &
787 e1_hat(2) = g2(3)*e3_hat(1) &
789 e1_hat(3) = g2(1)*e3_hat(2) &
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
798 e2_hat(1) = e3_hat(2)*e1_hat(3) &
800 e2_hat(2) = e3_hat(3)*e1_hat(1) &
802 e2_hat(3) = e3_hat(1)*e1_hat(2) &
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
814 (gausses(lx), shell, d, &
815 e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
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
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
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)
908 do jsize = 1, ndof*nn
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)
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)
927 xxi_lx = xi_lx /dsqrt( 1.0d0/3.0d0 )
928 eeta_lx = eta_lx/dsqrt( 3.0d0/5.0d0 )
930 do ip = 1, npoints_tying(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 ) )
941 xxi_lx = xi_lx /dsqrt( 3.0d0/5.0d0 )
942 eeta_lx = eta_lx/dsqrt( 1.0d0/3.0d0 )
944 do ip = 1, npoints_tying(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 ) )
955 xxi_lx = xi_lx /dsqrt( 1.0d0/3.0d0 )
956 eeta_lx = eta_lx/dsqrt( 1.0d0/3.0d0 )
958 do ip = 1, npoints_tying(3)
961 = ( 0.5d0*( 1.0d0+xxi_di(ip, 1)*xxi_lx ) ) &
962 *( 0.5d0*( 1.0d0+eeta_di(ip, 1)*eeta_lx ) )
966 do jsize = 1, ndof*nn
974 do ip = 1, npoints_tying(1)
978 = b(1, jsize)+h(ip, 1)*b_di(1, jsize, ip, 1, ly)
981 = b(5, jsize)+h(ip, 1)*b_di(5, jsize, ip, 1, ly)
985 do ip = 1, npoints_tying(2)
989 = b(2, jsize)+h(ip, 2)*b_di(2, jsize, ip, 2, ly)
992 = b(4, jsize)+h(ip, 2)*b_di(4, jsize, ip, 2, ly)
996 do ip = 1, npoints_tying(3)
1000 = b(3, jsize)+h(ip, 3)*b_di(3, jsize, ip, 3, ly)
1009 do jsize = 1, ndof*nn
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) )
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) )
1034 w_w_w_det = w_w_lx*w_ly*det
1038 db(1:5, 1:ndof*nn) = matmul( d, b(1:5, 1:ndof*nn ) )
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) )
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
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
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
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
1122 do jsize = 1, ndof*nn
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) )
1170 jsize = ndof*(nb-1)+j
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)
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
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)
1204 do jsize = 1, ndof*nn
1206 cv(jsize) = cv_theta(jsize)-0.5d0*cv_w(jsize)
1213 do jsize = 1, ndof*nn
1214 do isize = 1, ndof*nn
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)
1234 stiff(1:nn*ndof, 1:nn*ndof) = tmpstiff(1:nn*ndof, 1:nn*ndof)
1245 if( mixflag == 1 )
then
1273 tmpstiff(1:nn*ndof, 1:nn*ndof) = stiff(1:nn*ndof, 1:nn*ndof)
1277 stiff(i,j) = tmpstiff(sstable(i),sstable(j))
1281 elseif( mixflag == 2 )
then
1302 tmpstiff(1:nn*ndof, 1:nn*ndof) = stiff(1:nn*ndof, 1:nn*ndof)
1306 stiff(i,j) = tmpstiff(sstable(i),sstable(j))
1321 (etype, nn, ndof, ecoord, gausses, edisp, &
1322 strain, stress, thick, zeta, n_layer, n_totlyr)
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)
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
1348 integer :: npoints_tying(3)
1351 integer :: n_layer, n_totlyr
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), &
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), &
1388 real(kind = kreal) :: e(3, 3), ev(5)
1389 real(kind = kreal) :: s(3, 3), sv(5)
1390 real(kind = kreal) :: sumlyr
1429 elem(:, :) = ecoord(:, :)
1435 theta(1, na) = edisp(4, na)
1436 theta(2, na) = edisp(5, na)
1437 theta(3, na) = edisp(6, na)
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
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
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 )
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 )
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 )
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 )
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 )
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 )
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
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
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
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
1553 tpcoord(1, 1, 1) = 0.5d0
1554 tpcoord(2, 1, 1) = 0.0d0
1555 tpcoord(3, 1, 1) = 0.5d0
1557 tpcoord(1, 2, 1) = 0.0d0
1558 tpcoord(2, 2, 1) = 0.5d0
1559 tpcoord(3, 2, 1) = 0.5d0
1569 naturalcoord(1) = 0.0d0
1570 naturalcoord(2) = 0.0d0
1583 g1(i) = g1(i)+shapederiv(na, 1) &
1600 naturalcoord(1) = nncoord(nb, 1)
1601 naturalcoord(2) = nncoord(nb, 2)
1615 g1(i) = g1(i)+shapederiv(na, 1) &
1617 g2(i) = g2(i)+shapederiv(na, 2) &
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)
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) )
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
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)
1644 v2_abs = dsqrt( v2(1, nb)*v2(1, nb) &
1645 +v2(2, nb)*v2(2, nb) &
1646 +v2(3, nb)*v2(3, nb) )
1648 if( v2_abs .GT. 1.0d-15 )
then
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
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)
1661 v1_abs = dsqrt( v1(1, nb)*v1(1, nb) &
1662 +v1(2, nb)*v1(2, nb) &
1663 +v1(3, nb)*v1(3, nb) )
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
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)
1690 v3_abs = dsqrt( v3(1, nb)*v3(1, nb) &
1691 +v3(2, nb)*v3(2, nb) &
1692 +v3(3, nb)*v3(3, nb) )
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
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)
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)
1744 do ip = 1, npoints_tying(it)
1748 naturalcoord(1) = tpcoord(ip, 1, it)
1749 naturalcoord(2) = tpcoord(ip, 2, it)
1763 *( zeta_ly*a_over_2_v3(i, na) )
1766 = shapederiv(na, 1) &
1767 *( zeta_ly*a_over_2_v3(i, na) )
1769 = shapederiv(na, 2) &
1770 *( zeta_ly*a_over_2_v3(i, na) )
1771 dudzeta_rot(i, na) &
1773 *( a_over_2_v3(i, na) )
1790 g1(i) = g1(i)+shapederiv(na, 1) &
1793 g2(i) = g2(i)+shapederiv(na, 2) &
1796 g3(i) = g3(i)+dudzeta_rot(i, na)
1814 +shapederiv(na, 1) &
1816 +zeta_ly*a_over_2_theta_cross_v3(i, na) )
1819 +shapederiv(na, 2) &
1821 +zeta_ly*a_over_2_theta_cross_v3(i, na) )
1825 *( a_over_2_theta_cross_v3(i, na) )
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) ) )
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) ) )
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) )
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) )
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) )
1867 sumlyr = sumlyr +2*gausses(1)%pMaterial%shell_var(m)%weight
1869 zeta_ly = -1 +sumlyr -gausses(1)%pMaterial%shell_var(n_layer)%weight*(1-zeta)
1877 naturalcoord(1) = nncoord(lx, 1)
1878 naturalcoord(2) = nncoord(lx, 2)
1880 xi_lx = naturalcoord(1)
1881 eta_lx = naturalcoord(2)
1895 *( zeta_ly*a_over_2_v3(i, na) )
1898 = shapederiv(na, 1) &
1899 *( zeta_ly*a_over_2_v3(i, na) )
1901 = shapederiv(na, 2) &
1902 *( zeta_ly*a_over_2_v3(i, na) )
1903 dudzeta_rot(i, na) &
1905 *( a_over_2_v3(i, na) )
1922 g1(i) = g1(i)+shapederiv(na, 1) &
1925 g2(i) = g2(i)+shapederiv(na, 2) &
1928 g3(i) = g3(i)+dudzeta_rot(i, na)
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) )
1940 if(det == 0.0d0)
then
1941 write(*,*)
"ERROR:LIB Shell in l2009 Not Jacobian"
1951 *( g2(2)*g3(3)-g2(3)*g3(2) )
1953 *( g2(3)*g3(1)-g2(1)*g3(3) )
1955 *( g2(1)*g3(2)-g2(2)*g3(1) )
1957 *( g3(2)*g1(3)-g3(3)*g1(2) )
1959 *( g3(3)*g1(1)-g3(1)*g1(3) )
1961 *( g3(1)*g1(2)-g3(2)*g1(1) )
1963 *( g1(2)*g2(3)-g1(3)*g2(2) )
1965 *( g1(3)*g2(1)-g1(1)*g2(3) )
1967 *( g1(1)*g2(2)-g1(2)*g2(1) )
1971 g3_abs = dsqrt( g3(1)*g3(1) &
1979 e3_hat(1) = g3(1)/g3_abs
1980 e3_hat(2) = g3(2)/g3_abs
1981 e3_hat(3) = g3(3)/g3_abs
1983 e1_hat(1) = g2(2)*e3_hat(3) &
1985 e1_hat(2) = g2(3)*e3_hat(1) &
1987 e1_hat(3) = g2(1)*e3_hat(2) &
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
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
2021 +shapederiv(na, 1) &
2023 +zeta_ly*a_over_2_theta_cross_v3(i, na) )
2026 +shapederiv(na, 2) &
2028 +zeta_ly*a_over_2_theta_cross_v3(i, na) )
2032 *( a_over_2_theta_cross_v3(i, na) )
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) ) )
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) ) )
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) )
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) )
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) )
2072 = 0.5d0*( 1.0d0-xi_lx )*e23_di_2(4, 1) &
2073 +0.5d0*( 1.0d0+xi_lx )*e23_di_2(2, 1)
2076 = 0.5d0*( 1.0d0-eta_lx )*e31_di_2(1, 1) &
2077 +0.5d0*( 1.0d0+eta_lx )*e31_di_2(3, 1)
2082 xxi_lx = xi_lx /dsqrt( 1.0d0/3.0d0 )
2083 eeta_lx = eta_lx/dsqrt( 3.0d0/5.0d0 )
2085 do ip = 1, npoints_tying(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 ) )
2096 xxi_lx = xi_lx /dsqrt( 3.0d0/5.0d0 )
2097 eeta_lx = eta_lx/dsqrt( 1.0d0/3.0d0 )
2099 do ip = 1, npoints_tying(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 ) )
2110 xxi_lx = xi_lx /dsqrt( 1.0d0/3.0d0 )
2111 eeta_lx = eta_lx/dsqrt( 1.0d0/3.0d0 )
2113 do ip = 1, npoints_tying(3)
2116 = ( 0.5d0*( 1.0d0+xxi_di(ip, 1)*xxi_lx ) ) &
2117 *( 0.5d0*( 1.0d0+eeta_di(ip, 1)*eeta_lx ) )
2125 do ip = 1, npoints_tying(1)
2127 e11 = e11 +h(ip, 1)*e11_di(ip, 1)
2128 e31_2 = e31_2+h(ip, 1)*e31_di_2(ip, 1)
2136 do ip = 1, npoints_tying(2)
2138 e22 = e22 +h(ip, 2)*e22_di(ip, 2)
2139 e23_2 = e23_2+h(ip, 2)*e23_di_2(ip, 2)
2146 do ip = 1, npoints_tying(3)
2148 e12_2 = e12_2+h(ip, 3)*e12_di_2(ip, 3)
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) )
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) )
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)
2195 (gausses(lx), shell, d, &
2196 e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
2201 sv = matmul( d, ev )
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) )
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) )
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) )
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) )
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) )
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) )
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) )
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) )
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) )
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) )
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) )
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) )
2341 (etype, nn, ndof, xx, yy, zz, rho, thick, &
2342 ltype, params, vect, nsize, gausses)
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
2365 integer :: ivol, isurf
2372 integer :: jsize1, jsize2, jsize3, &
2373 jsize4, jsize5, jsize6
2375 integer :: n_totlyr, n_layer
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), &
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
2461 elem(1, na) = xx(na)
2462 elem(2, na) = yy(na)
2463 elem(3, na) = zz(na)
2476 do isize = 1, ndof*nn
2486 naturalcoord(1) = 0.0d0
2487 naturalcoord(2) = 0.0d0
2500 g1(i) = g1(i)+shapederiv(na, 1) &
2517 naturalcoord(1) = nncoord(nb, 1)
2518 naturalcoord(2) = nncoord(nb, 2)
2532 g1(i) = g1(i)+shapederiv(na, 1) &
2534 g2(i) = g2(i)+shapederiv(na, 2) &
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)
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) )
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
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)
2561 v2_abs = dsqrt( v2(1, nb)*v2(1, nb) &
2562 +v2(2, nb)*v2(2, nb) &
2563 +v2(3, nb)*v2(3, nb) )
2565 if( v2_abs .GT. 1.0d-15 )
then
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
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)
2578 v1_abs = dsqrt( v1(1, nb)*v1(1, nb) &
2579 +v1(2, nb)*v1(2, nb) &
2580 +v1(3, nb)*v1(3, nb) )
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
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)
2607 v3_abs = dsqrt( v3(1, nb)*v3(1, nb) &
2608 +v3(2, nb)*v3(2, nb) &
2609 +v3(3, nb)*v3(3, nb) )
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
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)
2632 if( ltype .LT. 10 )
then
2636 else if( ltype .GE. 10 )
then
2645 if( isurf .EQ. 1 )
then
2655 xi_lx = naturalcoord(1)
2656 eta_lx = naturalcoord(2)
2672 *( 0.0d0*a_over_2_v3(i, na) )
2675 = shapederiv(na, 1) &
2676 *( 0.0d0*a_over_2_v3(i, na) )
2678 = shapederiv(na, 2) &
2679 *( 0.0d0*a_over_2_v3(i, na) )
2680 dudzeta_rot(i, na) &
2682 *( a_over_2_v3(i, na) )
2699 g1(i) = g1(i)+shapederiv(na, 1) &
2702 g2(i) = g2(i)+shapederiv(na, 2) &
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)
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
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
2768 do isize = 1, ndof*nn
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
2789 if( ivol .EQ. 1 )
then
2792 n_totlyr = gausses(1)%pMaterial%totallyr
2793 do n_layer=1,n_totlyr
2801 sumlyr = sumlyr + 2 * gausses(1)%pMaterial%shell_var(m)%weight
2803 zeta_ly = -1 + sumlyr - gausses(1)%pMaterial%shell_var(n_layer)%weight*(1-
xg(ny, ly))
2816 xi_lx = naturalcoord(1)
2817 eta_lx = naturalcoord(2)
2833 *( zeta_ly*a_over_2_v3(i, na) )
2836 = shapederiv(na, 1) &
2837 *( zeta_ly*a_over_2_v3(i, na) )
2839 = shapederiv(na, 2) &
2840 *( zeta_ly*a_over_2_v3(i, na) )
2841 dudzeta_rot(i, na) &
2843 *( a_over_2_v3(i, na) )
2860 g1(i) = g1(i)+shapederiv(na, 1) &
2863 g2(i) = g2(i)+shapederiv(na, 2) &
2866 g3(i) = g3(i)+dudzeta_rot(i, na)
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) )
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
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
2914 w_w_w_det = w_w_lx*w_ly*det
2918 if( ltype .EQ. 1 )
then
2920 do isize = 1, ndof*nn
2922 vect(isize) = vect(isize)+w_w_w_det*n(1, isize)*val
2926 else if( ltype .EQ. 2 )
then
2928 do isize = 1, ndof*nn
2930 vect(isize) = vect(isize)+w_w_w_det*n(2, isize)*val
2934 else if( ltype .EQ. 3 )
then
2936 do isize = 1, ndof*nn
2938 vect(isize) = vect(isize)+w_w_w_det*n(3, isize)*val
2942 else if( ltype .EQ. 4 )
then
2944 do isize = 1, ndof*nn
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
2952 else if( ltype .EQ. 5 )
then
2960 x = x+shapefunc(nb)*elem(1, nb)
2961 y = y+shapefunc(nb)*elem(2, nb)
2962 z = z+shapefunc(nb)*elem(3, nb)
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
2974 coefx = phx*val*rho*val
2975 coefy = phy*val*rho*val
2976 coefz = phz*val*rho*val
2978 do isize = 1, ndof*nn
2982 +w_w_w_det*( n(1, isize)*coefx &
2983 +n(2, isize)*coefy &
2984 +n(3, isize)*coefz )
3016 (ic_type, nn, ndof, xx, yy, zz, rho, thick, &
3017 ltype, params, vect, nsize, gausses)
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
3038 real(kind = kreal) :: tmp(24)
3041 if(ic_type == 761)
then
3045 call dl_shell(731, 3, 6, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
3074 elseif(ic_type == 781)
then
3078 call dl_shell(741, 4, 6, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
3120 (etype, nn, ndof, ecoord, u, du, gausses, qf, thick, mixflag, nddisp)
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(:, :)
3136 real(kind = kreal),
intent(out) :: qf(:)
3137 real(kind = kreal),
intent(in) :: thick
3139 real(kind = kreal),
intent(in),
optional :: nddisp(3, nn)
3142 real(kind = kreal) :: stiff(nn*ndof, nn*ndof), totaldisp(nn*ndof)
3143 integer(kind = kint) :: i
3145 call stf_shell_mitc(etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag, nddisp)
3149 totaldisp(ndof*(i-1)+1:ndof*i) = u(1:ndof,i) + du(1:ndof,i)
3152 qf = matmul(stiff,totaldisp)
3159 (etype, nn, ndof, ecoord, u, du, gausses, qf, thick, mixflag, nddisp)
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)
3175 real(kind = kreal),
intent(out) :: qf(:)
3176 real(kind = kreal),
intent(in) :: thick
3178 real(kind = kreal),
intent(in),
optional :: nddisp(3, nn)
3181 real(kind = kreal) :: stiff(nn*ndof, nn*ndof), totaldisp(nn*ndof)
3182 integer(kind = kint) :: i
3184 call stf_shell_mitc(etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag, nddisp)
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)
3192 qf = matmul(stiff,totaldisp)
3197 subroutine mass_shell(etype, nn, elem, rho, thick, gausses, mass, lumped)
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
3209 real(kind=kreal),
intent(out) :: mass(:,:)
3210 real(kind=kreal),
intent(out) :: lumped(:)
3214 integer :: lx, ly, nsize, ndof
3219 integer :: jsize1, jsize2, jsize3, jsize4, jsize5, jsize6
3220 integer :: n_totlyr, n_layer
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
3272 naturalcoord(1) = 0.0d0
3273 naturalcoord(2) = 0.0d0
3280 g1(:) = matmul( elem, shapederiv(:,1) )
3290 naturalcoord(1) = nncoord(nb, 1)
3291 naturalcoord(2) = nncoord(nb, 2)
3297 g1(:) = matmul( elem, shapederiv(:,1) )
3298 g2(:) = matmul( elem, shapederiv(:,2) )
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)
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) )
3310 v3(:, nb) = det_cg3(:)/det_cg3_abs
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)
3318 v2_abs = dsqrt( v2(1, nb)*v2(1, nb) &
3319 +v2(2, nb)*v2(2, nb) &
3320 +v2(3, nb)*v2(3, nb) )
3322 if( v2_abs > 1.0d-15 )
then
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
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)
3335 v1_abs = dsqrt( v1(1, nb)*v1(1, nb) &
3336 +v1(2, nb)*v1(2, nb) &
3337 +v1(3, nb)*v1(3, nb) )
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
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)
3364 v3_abs = dsqrt( v3(1, nb)*v3(1, nb) &
3365 +v3(2, nb)*v3(2, nb) &
3366 +v3(3, nb)*v3(3, nb) )
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
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)
3386 n_totlyr = gausses(1)%pMaterial%totallyr
3387 do n_layer=1,n_totlyr
3394 sumlyr = sumlyr + 2 * gausses(1)%pMaterial%shell_var(m)%weight
3396 zeta_ly = -1 + sumlyr - gausses(1)%pMaterial%shell_var(n_layer)%weight*(1-
xg(ny, ly))
3409 xi_lx = naturalcoord(1)
3410 eta_lx = naturalcoord(2)
3423 u_rot(i, na) = shapefunc(na)*( zeta_ly*a_over_2_v3(i, na) )
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) )
3444 g1(i) = g1(i)+shapederiv(na, 1) *elem(i, na) &
3446 g2(i) = g2(i)+shapederiv(na, 2) *elem(i, na) &
3448 g3(i) = g3(i)+dudzeta_rot(i, na)
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) )
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
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
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
3508 totalmass = totalmass*3.d0
3514 totdiag = totdiag + mass(lx,lx)
3521 lumped(lx) = mass(lx,lx)/totdiag* totalmass
This module encapsulate the basic functions of all elements provide by this software.
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate deivatives of shape fucntion in natural coordiante system.
integer, parameter fe_mitc4_shell
integer, parameter fe_mitc9_shell
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
integer, parameter fe_mitc3_shell
subroutine getnodalnaturalcoord(fetype, nncoord)
This module provides data for gauss quadrature.
real(kind=kreal), dimension(3, 3) wgt
wieght of gauss points
real(kind=kreal), dimension(3, 3) xg
abscissa of gauss points
This module manages calculation relates with materials.
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.
This modules defines a structure to record history dependent parameter in static analysis.
All data should be recorded in every quadrature points.