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