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