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