FrontISTR  5.9.0
Large-scale structural analysis program with finit element method
static_LIB_shell.f90
Go to the documentation of this file.
1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2019 FrontISTR Commons
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
6 module m_static_lib_shell
7  use hecmw, only : kint, kreal
8  use elementinfo
9  use m_fstr_finiterotationkinematics, only: fstr_is_finite_rotation_shell_element
10 
11  !--------------------------------------------------------------------
12 
13  implicit none
14 
15  !--------------------------------------------------------------------
16 
17  !--------------------------------------------------------------
18  !
19  ! (Programmer)
20  ! Gaku Hashimoto
21  ! Department of Human and Engineered Environmental Studies
22  ! Graduate School of Frontier Sciences, The University of Tokyo
23  ! 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8563 JAPAN
24  !
25  ! (Ref.)
26  ! [1] Noguchi, H. and Hisada, T.,
27  ! "Sensitivity analysis in post-buckling problems of shell
28  ! structures,"
29  ! Computers & Structures, Vol.47, No.4, pp.699-710, (1993).
30  ! [2] Dvorkin, E.N. and Bathe, K.J.,
31  ! "A Continuum Mechanics Based Four-node Shell Element for
32  ! General Non-linear Analysis,"
33  ! Engineering Computations, Vol.1, pp.77-88, (1984).
34  ! [3] Bucalem, M.L. and Bathe, K.J.,
35  ! "Higher-order MITC general shell element,"
36  ! International Journal for Numerical Methods in
37  ! Engineering, Vol.36, pp.3729-3754, (1993).
38  ! [4] Lee, P.S. and Bathe, K.J.,
39  ! "Development of MITC Isotropic Triangular Shell Finite
40  ! Elements,"
41  ! Computers & Structures, Vol.82, pp.945-962, (2004).
42  !
43  ! Xi YUAN
44  ! Apr. 13, 2019: Introduce mass matrix calculation
45  ! (Ref.)
46  ! [5] E.Hinton, T.A.Rock, O.C.Zienkiewicz(1976): A Note on Mass Lumping
47  ! and Related Process in FEM. International Journal on Earthquake Eng
48  ! and structural dynamics, 4, pp245-249
49  !
50  !--------------------------------------------------------------
51 
52  !--------------------------------------------------------------------
53 
54 contains
55 
56  logical function shellsupportsfiniterotationkinematics( etype, nn )
57  implicit none
58  integer(kind=kint), intent(in) :: etype, nn
59 
60  shellsupportsfiniterotationkinematics = fstr_is_finite_rotation_shell_element( etype, nn )
61  end function shellsupportsfiniterotationkinematics
62 
63  pure function outer_product3(a, b) result(ab)
64 
65  real(kind = kreal), intent(in) :: a(3), b(3)
66  real(kind = kreal) :: ab(3, 3)
67  integer :: i, j
68 
69  do j = 1, 3
70  do i = 1, 3
71  ab(i, j) = a(i)*b(j)
72  end do
73  end do
74 
75  end function outer_product3
76 
77 
78  pure subroutine shell_basis_from_covariant(g1, g2, g3, e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, det)
79 
80  real(kind = kreal), intent(in) :: g1(3), g2(3), g3(3)
81  real(kind = kreal), intent(out) :: e1_hat(3), e2_hat(3), e3_hat(3)
82  real(kind = kreal), intent(out) :: cg1(3), cg2(3), cg3(3)
83  real(kind = kreal), intent(out) :: det
84  real(kind = kreal) :: det_inv, g3_abs, e1_abs, e2_abs
85 
86  det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
87  +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
88  +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )
89  det_inv = 1.0d0/det
90 
91  cg1(1) = det_inv*( g2(2)*g3(3)-g2(3)*g3(2) )
92  cg1(2) = det_inv*( g2(3)*g3(1)-g2(1)*g3(3) )
93  cg1(3) = det_inv*( g2(1)*g3(2)-g2(2)*g3(1) )
94  cg2(1) = det_inv*( g3(2)*g1(3)-g3(3)*g1(2) )
95  cg2(2) = det_inv*( g3(3)*g1(1)-g3(1)*g1(3) )
96  cg2(3) = det_inv*( g3(1)*g1(2)-g3(2)*g1(1) )
97  cg3(1) = det_inv*( g1(2)*g2(3)-g1(3)*g2(2) )
98  cg3(2) = det_inv*( g1(3)*g2(1)-g1(1)*g2(3) )
99  cg3(3) = det_inv*( g1(1)*g2(2)-g1(2)*g2(1) )
100 
101  g3_abs = dsqrt( dot_product(g3, g3) )
102  e3_hat(:) = g3(:)/g3_abs
103  e1_hat(1) = g2(2)*e3_hat(3)-g2(3)*e3_hat(2)
104  e1_hat(2) = g2(3)*e3_hat(1)-g2(1)*e3_hat(3)
105  e1_hat(3) = g2(1)*e3_hat(2)-g2(2)*e3_hat(1)
106  e1_abs = dsqrt( dot_product(e1_hat, e1_hat) )
107  e1_hat(:) = e1_hat(:)/e1_abs
108  e2_hat(1) = e3_hat(2)*e1_hat(3)-e3_hat(3)*e1_hat(2)
109  e2_hat(2) = e3_hat(3)*e1_hat(1)-e3_hat(1)*e1_hat(3)
110  e2_hat(3) = e3_hat(1)*e1_hat(2)-e3_hat(2)*e1_hat(1)
111  e2_abs = dsqrt( dot_product(e2_hat, e2_hat) )
112  e2_hat(:) = e2_hat(:)/e2_abs
113 
114  end subroutine shell_basis_from_covariant
115 
116  !--------------------------------------------------------------------
117  subroutine shellstressvectortotensor( stress, tensor )
118  implicit none
119 
120  real(kind=kreal), intent(in) :: stress(6)
121  real(kind=kreal), intent(out) :: tensor(3, 3)
122 
123  tensor(:, :) = 0.0d0
124  tensor(1, 1) = stress(1)
125  tensor(2, 2) = stress(2)
126  tensor(3, 3) = stress(3)
127  tensor(1, 2) = stress(4)
128  tensor(2, 1) = tensor(1, 2)
129  tensor(2, 3) = stress(5)
130  tensor(3, 2) = tensor(2, 3)
131  tensor(3, 1) = stress(6)
132  tensor(1, 3) = tensor(3, 1)
133 
134  end subroutine shellstressvectortotensor
135 
136  subroutine shelltensortostressvector( tensor, stress )
137  implicit none
138 
139  real(kind=kreal), intent(in) :: tensor(3, 3)
140  real(kind=kreal), intent(out) :: stress(6)
141 
142  stress(1) = tensor(1, 1)
143  stress(2) = tensor(2, 2)
144  stress(3) = tensor(3, 3)
145  stress(4) = tensor(1, 2)
146  stress(5) = tensor(2, 3)
147  stress(6) = tensor(3, 1)
148 
149  end subroutine shelltensortostressvector
150 
151  !--------------------------------------------------------------------
152  subroutine shellobjectivestressincrement( stress_old, dstrain, dstress_obj, trace_coeff )
153  implicit none
154 
155  real(kind=kreal), intent(in) :: stress_old(6), dstrain(6)
156  real(kind=kreal), intent(out) :: dstress_obj(6)
157  real(kind=kreal), intent(in), optional :: trace_coeff
158 
159  call shellobjectivetracestressincrement( stress_old, dstrain, dstress_obj, trace_coeff )
160 
161  end subroutine shellobjectivestressincrement
162 
163  !--------------------------------------------------------------------
164  subroutine shellobjectivetracestressincrement( stress_old, dstrain, dstress_trace, trace_coeff )
165  implicit none
166 
167  real(kind=kreal), intent(in) :: stress_old(6), dstrain(6)
168  real(kind=kreal), intent(out) :: dstress_trace(6)
169  real(kind=kreal), intent(in), optional :: trace_coeff
170 
171  real(kind=kreal) :: stress_tensor(3, 3), dstress_tensor(3, 3)
172  real(kind=kreal) :: trd, coeff
173 
174  call shellstressvectortotensor( stress_old, stress_tensor )
175  coeff = 1.0d0
176  if( present( trace_coeff ) ) coeff = trace_coeff
177  trd = coeff*(dstrain(1)+dstrain(2))+dstrain(3)
178  dstress_tensor(:, :) = -stress_tensor(:, :)*trd
179  call shelltensortostressvector( dstress_tensor, dstress_trace )
180 
181  end subroutine shellobjectivetracestressincrement
182 
183  !--------------------------------------------------------------------
184  real(kind=kreal) function shellplanestresstracecoeff( gauss, n_layer )
185  use mmechgauss
186  use mmaterial, only: m_poisson, mc_isoelastic, getelastictype
187  implicit none
188 
189  type(tgaussstatus), intent(in) :: gauss
190  integer(kind=kint), intent(in) :: n_layer
191 
192  real(kind=kreal) :: nu, outa(2)
193  logical :: ierr
194 
195  ! Isotropic plane-stress trace coefficient for the UL stress update.
196  shellplanestresstracecoeff = 1.0d0
197  if( .not. associated( gauss%pMaterial ) ) return
198  if( getelastictype( gauss%pMaterial%mtype ) == 1 ) then
199  stop "MITC4 shell UL orthotropic trace correction is not supported"
200  endif
201  nu = gauss%pMaterial%variables(m_poisson)
202  call fetch_tabledata(mc_isoelastic, gauss%pMaterial%dict, outa, ierr)
203  if( associated( gauss%pMaterial%shell_var ) ) then
204  if( n_layer >= 1 .and. n_layer <= size( gauss%pMaterial%shell_var ) ) then
205  if( gauss%pMaterial%shell_var(n_layer)%ortho == 0 ) then
206  if( ierr ) then
207  nu = gauss%pMaterial%shell_var(n_layer)%pp
208  else
209  nu = outa(2)
210  endif
211  else
212  stop "MITC4 shell UL orthotropic trace correction is not supported"
213  endif
214  else if( .not. ierr ) then
215  nu = outa(2)
216  endif
217  else if( .not. ierr ) then
218  nu = outa(2)
219  endif
220 
221  if( abs(1.0d0-nu) > 1.0d-12 ) then
222  shellplanestresstracecoeff = (1.0d0-2.0d0*nu)/(1.0d0-nu)
223  endif
224 
225  end function shellplanestresstracecoeff
226 
227  !--------------------------------------------------------------------
228  subroutine shelladdstressvectortodb( dstress, j, DB )
229  implicit none
230 
231  integer(kind=kint), intent(in) :: j
232  real(kind=kreal), intent(in) :: dstress(6)
233  real(kind=kreal), intent(inout) :: db(:, :)
234 
235  db(1, j) = db(1, j)+dstress(1)
236  db(2, j) = db(2, j)+dstress(2)
237  db(3, j) = db(3, j)+dstress(4)
238  db(4, j) = db(4, j)+dstress(5)
239  db(5, j) = db(5, j)+dstress(6)
240 
241  end subroutine shelladdstressvectortodb
242 
243  !--------------------------------------------------------------------
244  subroutine shelladdulobjectivetracetangent( stress_old, ncol, B, DB, trace_coeff )
245  implicit none
246 
247  integer(kind=kint), intent(in) :: ncol
248  real(kind=kreal), intent(in) :: stress_old(6)
249  real(kind=kreal), intent(in) :: b(5, ncol)
250  real(kind=kreal), intent(inout) :: db(5, ncol)
251  real(kind=kreal), intent(in), optional :: trace_coeff
252 
253  integer(kind=kint) :: j
254  real(kind=kreal) :: dstrain_col(6), dstress_trace(6)
255 
256  do j = 1, ncol
257  dstrain_col(:) = 0.0d0
258  dstrain_col(1) = b(1, j)
259  dstrain_col(2) = b(2, j)
260  call shellobjectivetracestressincrement( stress_old, dstrain_col, dstress_trace, trace_coeff )
261  call shelladdstressvectortodb( dstress_trace, j, db )
262  end do
263 
264  end subroutine shelladdulobjectivetracetangent
265 
266  !--------------------------------------------------------------------
267  subroutine shelladdulobjectivetangent( stress_old, ncol, B, DB, trace_coeff )
268  implicit none
269 
270  integer(kind=kint), intent(in) :: ncol
271  real(kind=kreal), intent(in) :: stress_old(6)
272  real(kind=kreal), intent(in) :: b(5, ncol)
273  real(kind=kreal), intent(inout) :: db(5, ncol)
274  real(kind=kreal), intent(in), optional :: trace_coeff
275 
276  call shelladdulobjectivetracetangent( stress_old, ncol, b, db, trace_coeff )
277 
278  end subroutine shelladdulobjectivetangent
279 
280 
281  !####################################################################
282  subroutine stf_shell_mitc &
283  (etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag, nddisp, element, qf_stress, include_geo_stiff, &
284  nddirector, ndrefdirector, ndcurdirector, nddrill)
285  !####################################################################
286 
287  use mmechgauss
288  use quadrature
289  use m_matmatrix
290  use mmaterial, only: infinitesimal, updatelag, totallag, iselastic
291 
292  !--------------------------------------------------------------------
293 
294  integer(kind = kint), intent(in) :: etype
295  integer(kind = kint), intent(in) :: nn, mixflag
296  integer(kind = kint), intent(in) :: ndof
297  real(kind = kreal), intent(in) :: ecoord(3, nn)
298  type(tgaussstatus), intent(in) :: gausses(:)
299  real(kind = kreal), intent(out) :: stiff(:, :)
300  real(kind = kreal), intent(in) :: thick
301 
302  real(kind = kreal), intent(in), optional :: nddisp(ndof, nn)
303  type(telement), intent(in), optional :: element
304  real(kind = kreal), intent(out), optional :: qf_stress(:)
305  logical, intent(in), optional :: include_geo_stiff
306  real(kind = kreal), intent(in), optional :: nddirector(3, nn)
307  real(kind = kreal), intent(in), optional :: ndrefdirector(3, nn)
308  real(kind = kreal), intent(in), optional :: ndcurdirector(3, nn)
309  real(kind = kreal), intent(in), optional :: nddrill(nn)
310 
311  !--------------------------------------------------------------------
312 
313  integer :: flag, flag_dof
314  integer :: ndof_shell
315  integer(kind=kint) :: ierr_quad
316  integer(kind=kint) :: ishell
317  integer :: i, j, m, n
318  integer :: lx, ly
319  integer :: fetype
320  integer :: ny
321  integer :: ntying
322  integer :: npoints_tying(3)
323  integer :: it, ip
324  integer :: na, nb
325  integer :: isize, jsize
326  integer :: jsize1, jsize2, jsize3, &
327  jsize4, jsize5, jsize6
328  integer :: n_layer,n_totlyr, sstable(24)
329 
330  real(kind = kreal) :: d(5, 5), b(5, ndof*nn), db(5, ndof*nn)
331  real(kind = kreal) :: stress_old_vec(6)
332  real(kind = kreal) :: tmpstiff(ndof*nn, ndof*nn)
333  real(kind = kreal) :: qf_tmp(ndof*nn), qf_mix(ndof*nn), qf_disp(ndof*nn)
334  real(kind = kreal) :: qf_stress_integrand(ndof*nn), vol_deriv(ndof*nn)
335  real(kind = kreal) :: sv_force(5), geo_term
336  real(kind = kreal) :: s_global(3, 3), smat(9, 9)
337  real(kind = kreal) :: bn(9, ndof*nn), sbn(9, ndof*nn)
338  real(kind = kreal) :: elem(3, nn)
339  real(kind = kreal) :: shell_disp(6, nn)
340  real(kind = kreal) :: xi_lx, eta_lx, zeta_ly
341  real(kind = kreal) :: w_w_lx, w_ly
342  real(kind = kreal) :: b_di(5, ndof*nn, 6, 3, 7)
343  real(kind = kreal) :: bg1_di(3, ndof*nn, 6, 3, 7)
344  real(kind = kreal) :: bg2_di(3, ndof*nn, 6, 3, 7)
345  real(kind = kreal) :: bg3_di(3, ndof*nn, 6, 3, 7)
346  real(kind = kreal) :: b1(3, ndof*nn), b2(3, ndof*nn), &
347  b3(3, ndof*nn)
348  real(kind = kreal), allocatable :: b2rot(:, :, :, :)
349  real(kind = kreal), allocatable :: b2rot_di(:, :, :, :, :, :, :)
350  real(kind = kreal) :: naturalcoord(2)
351  real(kind = kreal) :: tpcoord(6, 2, 3)
352  real(kind = kreal) :: nncoord(nn, 2)
353  real(kind = kreal) :: shapefunc(nn)
354  real(kind = kreal) :: shapederiv(nn, 2)
355  real(kind = kreal) :: aa1(3), aa2(3), aa3(3)
356  real(kind = kreal) :: bb1(3), bb2(3), bb3(3)
357  real(kind = kreal) :: cc1(3), cc2(3)
358  real(kind = kreal) :: alpha
359  real(kind = kreal) :: trace_coeff
360  real(kind = kreal) :: xxi_lx, eeta_lx
361  real(kind = kreal) :: xxi_di(6, 3), eeta_di(6, 3)
362  real(kind = kreal) :: h(nn, 3)
363  real(kind = kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
364  real(kind = kreal) :: v1_i(3), v2_i(3), v3_i(3)
365  real(kind = kreal) :: v1_abs, v2_abs, v3_abs
366  real(kind = kreal) :: a_over_2_v3(3, nn)
367  real(kind = kreal) :: a_over_2_v3_ref(3, nn)
368  real(kind = kreal) :: a_over_2_v3_deriv(3, 3, nn)
369  real(kind = kreal) :: a_over_2_v3_second(3, 3, 3, nn)
370  real(kind = kreal) :: u_rot(3, nn)
371  real(kind = kreal) :: dudxi_rot(3, nn), dudeta_rot(3, nn), &
372  dudzeta_rot(3, nn)
373  real(kind = kreal) :: dudxi_rot_deriv(3, 3), dudeta_rot_deriv(3, 3), &
374  dudzeta_rot_deriv(3, 3)
375  real(kind = kreal) :: dudxi_rot_second(3, 3, 3), dudeta_rot_second(3, 3, 3), &
376  dudzeta_rot_second(3, 3, 3)
377  real(kind = kreal) :: g1(3), g2(3), g3(3)
378  real(kind = kreal) :: g1_tl(3), g2_tl(3)
379  real(kind = kreal) :: g1_weight(3), g2_weight(3), g3_weight(3)
380  real(kind = kreal) :: dudxi_trans(3), dudeta_trans(3)
381  real(kind = kreal) :: g3_abs
382  real(kind = kreal) :: e_0(3)
383  real(kind = kreal) :: cg1(3), cg2(3), cg3(3)
384  real(kind = kreal) :: det
385  real(kind = kreal) :: det_weight
386  real(kind = kreal) :: det_cg3(3)
387  real(kind = kreal) :: det_inv
388  real(kind = kreal) :: det_cg3_abs
389  real(kind = kreal) :: w_w_w_det
390  real(kind = kreal) :: e1_hat(3), e2_hat(3), e3_hat(3)
391  real(kind = kreal) :: e1_hat_abs, e2_hat_abs
392  real(kind = kreal) :: e1_hat_mat(3), e2_hat_mat(3), e3_hat_mat(3)
393  real(kind = kreal) :: cg1_mat(3), cg2_mat(3), cg3_mat(3)
394  real(kind = kreal) :: det_mat
395  real(kind = kreal) :: cv12(ndof*nn), cv13(ndof*nn), &
396  cv21(ndof*nn), cv23(ndof*nn), &
397  cv31(ndof*nn), cv32(ndof*nn)
398  real(kind = kreal) :: cv_theta(ndof*nn), cv_w(ndof*nn)
399  real(kind = kreal) :: cv(ndof*nn)
400  real(kind = kreal) :: cv_w_second(3, 3, nn)
401  real(kind = kreal) :: cv_disp, cv_deriv, cv_deriv_disp
402  real(kind = kreal) :: cv12_2, cv13_2, cv21_2, cv23_2, cv31_2, cv32_2
403  real(kind = kreal) :: drill_axis(3), rot_projector(3, 3), hrot(3, 3)
404  real(kind = kreal) :: hess_coeff, drill_coeff, axis_norm
405  real(kind = kreal) :: director_inc(3), director_ref(3), director_deriv(3, 3)
406  real(kind = kreal) :: director_second(3, 3, 3)
407  logical :: finite_rotation_director, add_geo_stiff, use_tl_green, use_director_tangent
408 
409  sstable = 0
410  flag_dof = 0
411  ny = 0
412  shell_disp(:, :) = 0.0d0
413  b_di(:, :, :, :, :) = 0.0d0
414  b1(:, :) = 0.0d0
415  b2(:, :) = 0.0d0
416  b3(:, :) = 0.0d0
417  a_over_2_v3_second(:, :, :, :) = 0.0d0
418  ndof_shell = min(ndof, 6)
419 
420  !--------------------------------------------------------------------
421 
422  ! MITC4
423  if( etype .EQ. fe_mitc4_shell ) then
424 
425  fetype = fe_mitc4_shell
426 
427  ny = 2
428 
429  ntying = 1
430  npoints_tying(1)= 4
431 
432  ! MITC9
433  else if( etype .EQ. fe_mitc9_shell ) then
434 
435  fetype = fe_mitc9_shell
436 
437  ny = 3
438 
439  ntying = 3
440  npoints_tying(1)= 6
441  npoints_tying(2)= 6
442  npoints_tying(3)= 4
443 
444  ! MITC3
445  else if( etype .EQ. fe_mitc3_shell ) then
446 
447  fetype = fe_mitc3_shell
448 
449  ny = 2
450 
451  ntying = 1
452  npoints_tying(1)= 3
453 
454  end if
455 
456  !--------------------------------------------------------------------
457 
458  if( present( nddisp ) ) then
459 
460  shell_disp(1:ndof_shell, 1:nn) = nddisp(1:ndof_shell, 1:nn)
461 
462  end if
463 
464  !--------------------------------------------------------------------
465 
466  flag = gausses(1)%pMaterial%nlgeom_flag
467 
468  if( .not. present( nddisp ) ) flag = infinitesimal
469  finite_rotation_director = ( flag == totallag .or. flag == updatelag ) .and. ndof_shell >= 6 &
470  .and. shellsupportsfiniterotationkinematics( etype, nn ) &
471  .and. iselastic(gausses(1)%pMaterial%mtype)
472  use_director_tangent = finite_rotation_director .and. iselastic(gausses(1)%pMaterial%mtype)
473  add_geo_stiff = flag /= infinitesimal
474  if( present( include_geo_stiff ) ) add_geo_stiff = add_geo_stiff .and. include_geo_stiff
475  ! Green-Lagrange strain for elastic finite-rotation shell elements.
476  use_tl_green = flag == totallag .and. shellsupportsfiniterotationkinematics( etype, nn ) &
477  .and. iselastic(gausses(1)%pMaterial%mtype)
478  if( use_director_tangent ) then
479  allocate( b2rot(5, 3, 3, nn) )
480  allocate( b2rot_di(5, 3, 3, nn, 6, 3, ny) )
481  b2rot_di(:, :, :, :, :, :, :) = 0.0d0
482  endif
483 
484  !--------------------------------------------------------------------
485 
486  elem(:, :) = ecoord(:, :)
487  if( flag == updatelag ) elem(:, :) = elem(:, :) + shell_disp(1:3, :)
488  bg1_di(:, :, :, :, :) = 0.0d0
489  bg2_di(:, :, :, :, :) = 0.0d0
490  bg3_di(:, :, :, :, :) = 0.0d0
491 
492  !--------------------------------------------------------------------
493 
494  tmpstiff(:, :) = 0.0d0
495  qf_tmp(:) = 0.0d0
496  qf_disp(:) = 0.0d0
497  do nb = 1, nn
498  qf_disp(ndof*(nb-1)+1:ndof*(nb-1)+ndof_shell) = shell_disp(1:ndof_shell, nb)
499  end do
500  if( present( qf_stress ) ) then
501  qf_stress(1:ndof*nn) = 0.0d0
502  endif
503 
504  !-------------------------------------------------------------------
505 
506  ! xi-coordinate at a node in a local element
507  ! eta-coordinate at a node in a local element
508  call getnodalnaturalcoord(fetype, nncoord)
509 
510  !-------------------------------------------------------------------
511 
512  ! MITC4
513  if( etype .EQ. fe_mitc4_shell ) then
514 
515  !--------------------------------------------------------
516 
517  ! xi-coordinate at a tying point in a local element
518  tpcoord(1, 1, 1) = 0.0d0
519  tpcoord(2, 1, 1) = 1.0d0
520  tpcoord(3, 1, 1) = 0.0d0
521  tpcoord(4, 1, 1) = -1.0d0
522  ! eta-coordinate at a tying point in a local element
523  tpcoord(1, 2, 1) = -1.0d0
524  tpcoord(2, 2, 1) = 0.0d0
525  tpcoord(3, 2, 1) = 1.0d0
526  tpcoord(4, 2, 1) = 0.0d0
527 
528  !--------------------------------------------------------
529 
530  ! MITC9
531  else if( etype .EQ. fe_mitc9_shell ) then
532 
533  !--------------------------------------------------------
534 
535  ! xi-coordinate at a tying point in a local element
536  tpcoord(1, 1, 1) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
537  tpcoord(2, 1, 1) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
538  tpcoord(3, 1, 1) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
539  tpcoord(4, 1, 1) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
540  tpcoord(5, 1, 1) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
541  tpcoord(6, 1, 1) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
542  ! eta-coordinate at a tying point in a local element
543  tpcoord(1, 2, 1) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
544  tpcoord(2, 2, 1) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
545  tpcoord(3, 2, 1) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
546  tpcoord(4, 2, 1) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
547  tpcoord(5, 2, 1) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
548  tpcoord(6, 2, 1) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
549 
550  ! xi-coordinate at a tying point in a local element
551  tpcoord(1, 1, 2) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
552  tpcoord(2, 1, 2) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
553  tpcoord(3, 1, 2) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
554  tpcoord(4, 1, 2) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
555  tpcoord(5, 1, 2) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
556  tpcoord(6, 1, 2) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
557  ! eta-coordinate at a tying point in a local element
558  tpcoord(1, 2, 2) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
559  tpcoord(2, 2, 2) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
560  tpcoord(3, 2, 2) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
561  tpcoord(4, 2, 2) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
562  tpcoord(5, 2, 2) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
563  tpcoord(6, 2, 2) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
564 
565  ! xi-coordinate at a tying point in a local element
566  tpcoord(1, 1, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
567  tpcoord(2, 1, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
568  tpcoord(3, 1, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
569  tpcoord(4, 1, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
570  ! eta-coordinate at a tying point in a local element
571  tpcoord(1, 2, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
572  tpcoord(2, 2, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
573  tpcoord(3, 2, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
574  tpcoord(4, 2, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
575 
576  !--------------------------------------------------------
577 
578  ! Xi-coordinate at a tying point in a local element
579  xxi_di(1, 1) = -1.0d0
580  xxi_di(2, 1) = 1.0d0
581  xxi_di(3, 1) = 1.0d0
582  xxi_di(4, 1) = -1.0d0
583  xxi_di(5, 1) = 1.0d0
584  xxi_di(6, 1) = -1.0d0
585  ! Eta-coordinate at a tying point in a local element
586  eeta_di(1, 1) = -1.0d0
587  eeta_di(2, 1) = -1.0d0
588  eeta_di(3, 1) = 1.0d0
589  eeta_di(4, 1) = 1.0d0
590  eeta_di(5, 1) = 0.0d0
591  eeta_di(6, 1) = 0.0d0
592 
593  ! Xi-coordinate at a tying point in a local element
594  xxi_di(1, 2) = -1.0d0
595  xxi_di(2, 2) = 0.0d0
596  xxi_di(3, 2) = 1.0d0
597  xxi_di(4, 2) = 1.0d0
598  xxi_di(5, 2) = 0.0d0
599  xxi_di(6, 2) = -1.0d0
600  ! Eta-coordinate at a tying point in a local element
601  eeta_di(1, 2) = -1.0d0
602  eeta_di(2, 2) = -1.0d0
603  eeta_di(3, 2) = -1.0d0
604  eeta_di(4, 2) = 1.0d0
605  eeta_di(5, 2) = 1.0d0
606  eeta_di(6, 2) = 1.0d0
607 
608  !--------------------------------------------------------
609 
610  ! MITC3
611  else if( etype .EQ. fe_mitc3_shell ) then
612 
613  !--------------------------------------------------------
614 
615  ! xi-coordinate at a tying point in a local element
616  tpcoord(1, 1, 1) = 0.5d0
617  tpcoord(2, 1, 1) = 0.0d0
618  tpcoord(3, 1, 1) = 0.5d0
619  ! eta-coordinate at a tying point in a local element
620  tpcoord(1, 2, 1) = 0.0d0
621  tpcoord(2, 2, 1) = 0.5d0
622  tpcoord(3, 2, 1) = 0.5d0
623 
624  !--------------------------------------------------------
625 
626  end if
627 
628  !--------------------------------------------------------------------
629 
630  ! xi-coordinate at the center point in a local element
631  ! eta-coordinate at the center point in a local element
632  naturalcoord(1) = 0.0d0
633  naturalcoord(2) = 0.0d0
634 
635  call getshapederiv(fetype, naturalcoord, shapederiv)
636 
637  !--------------------------------------------------------------
638 
639  ! Covariant basis vector
640  do i = 1, 3
641 
642  g1(i) = 0.0d0
643 
644  do na = 1, nn
645 
646  g1(i) = g1(i)+shapederiv(na, 1) &
647  *elem(i, na)
648 
649  end do
650 
651  end do
652 
653  e_0(1) = g1(1)
654  e_0(2) = g1(2)
655  e_0(3) = g1(3)
656 
657  !--------------------------------------------------------------
658 
659  do nb = 1, nn
660 
661  !--------------------------------------------------------
662 
663  naturalcoord(1) = nncoord(nb, 1)
664  naturalcoord(2) = nncoord(nb, 2)
665 
666  call getshapederiv(fetype, naturalcoord, shapederiv)
667 
668  !--------------------------------------------------------
669 
670  ! Covariant basis vector
671  do i = 1, 3
672 
673  g1(i) = 0.0d0
674  g2(i) = 0.0d0
675 
676  do na = 1, nn
677 
678  g1(i) = g1(i)+shapederiv(na, 1) &
679  *elem(i, na)
680  g2(i) = g2(i)+shapederiv(na, 2) &
681  *elem(i, na)
682 
683  end do
684 
685  end do
686 
687  !------------------------------------------
688 
689  det_cg3(1) = g1(2)*g2(3)-g1(3)*g2(2)
690  det_cg3(2) = g1(3)*g2(1)-g1(1)*g2(3)
691  det_cg3(3) = g1(1)*g2(2)-g1(2)*g2(1)
692 
693  det_cg3_abs = dsqrt( det_cg3(1)*det_cg3(1) &
694  +det_cg3(2)*det_cg3(2) &
695  +det_cg3(3)*det_cg3(3) )
696 
697  v3(1, nb) = det_cg3(1)/det_cg3_abs
698  v3(2, nb) = det_cg3(2)/det_cg3_abs
699  v3(3, nb) = det_cg3(3)/det_cg3_abs
700 
701  !--------------------------------------------------------
702 
703  v2(1, nb) = v3(2, nb)*e_0(3)-v3(3, nb)*e_0(2)
704  v2(2, nb) = v3(3, nb)*e_0(1)-v3(1, nb)*e_0(3)
705  v2(3, nb) = v3(1, nb)*e_0(2)-v3(2, nb)*e_0(1)
706 
707  v2_abs = dsqrt( v2(1, nb)*v2(1, nb) &
708  +v2(2, nb)*v2(2, nb) &
709  +v2(3, nb)*v2(3, nb) )
710 
711  if( v2_abs .GT. 1.0d-15 ) then
712 
713  v2(1, nb) = v2(1, nb)/v2_abs
714  v2(2, nb) = v2(2, nb)/v2_abs
715  v2(3, nb) = v2(3, nb)/v2_abs
716 
717  v1(1, nb) = v2(2, nb)*v3(3, nb) &
718  -v2(3, nb)*v3(2, nb)
719  v1(2, nb) = v2(3, nb)*v3(1, nb) &
720  -v2(1, nb)*v3(3, nb)
721  v1(3, nb) = v2(1, nb)*v3(2, nb) &
722  -v2(2, nb)*v3(1, nb)
723 
724  v1_abs = dsqrt( v1(1, nb)*v1(1, nb) &
725  +v1(2, nb)*v1(2, nb) &
726  +v1(3, nb)*v1(3, nb) )
727 
728  v1(1, nb) = v1(1, nb)/v1_abs
729  v1(2, nb) = v1(2, nb)/v1_abs
730  v1(3, nb) = v1(3, nb)/v1_abs
731 
732  else ! YX: impossible
733 
734  v1(1, nb) = 0.0d0
735  v1(2, nb) = 0.0d0
736  v1(3, nb) = -1.0d0
737 
738  v2(1, nb) = 0.0d0
739  v2(2, nb) = 1.0d0
740  v2(3, nb) = 0.0d0
741 
742  end if
743 
744  !---------------------------------------------------
745 
746  v3(1, nb) = v1(2, nb)*v2(3, nb) &
747  -v1(3, nb)*v2(2, nb)
748  v3(2, nb) = v1(3, nb)*v2(1, nb) &
749  -v1(1, nb)*v2(3, nb)
750  v3(3, nb) = v1(1, nb)*v2(2, nb) &
751  -v1(2, nb)*v2(1, nb)
752 
753  v3_abs = dsqrt( v3(1, nb)*v3(1, nb) &
754  +v3(2, nb)*v3(2, nb) &
755  +v3(3, nb)*v3(3, nb) )
756 
757  v3(1, nb) = v3(1, nb)/v3_abs
758  v3(2, nb) = v3(2, nb)/v3_abs
759  v3(3, nb) = v3(3, nb)/v3_abs
760 
761  !--------------------------------------------------------
762 
763  a_over_2_v3(1, nb) = 0.5d0*thick*v3(1, nb)
764  a_over_2_v3(2, nb) = 0.5d0*thick*v3(2, nb)
765  a_over_2_v3(3, nb) = 0.5d0*thick*v3(3, nb)
766  a_over_2_v3_ref(1:3, nb) = a_over_2_v3(1:3, nb)
767  if( finite_rotation_director .and. flag == updatelag .and. present( ndcurdirector ) ) then
768  a_over_2_v3_ref(1:3, nb) = ndcurdirector(1:3, nb)
769  else if( finite_rotation_director .and. present( ndrefdirector ) ) then
770  a_over_2_v3_ref(1:3, nb) = ndrefdirector(1:3, nb)
771  endif
772 
773  director_ref(1:3) = a_over_2_v3_ref(1:3, nb)
774  if( finite_rotation_director ) then
775  if( present( nddirector ) ) then
776  a_over_2_v3(1:3, nb) = nddirector(1:3, nb)
777  else
778  call shelldirectorincrement( shell_disp(4:6, nb), director_ref, director_inc )
779  a_over_2_v3(1:3, nb) = director_ref(1:3) + director_inc(1:3)
780  endif
781  if( use_director_tangent ) then
782  call shelldirectorincrementalderiv( a_over_2_v3(1:3, nb), director_deriv )
783  call shelldirectorincrementalsecondderiv( a_over_2_v3(1:3, nb), director_second )
784  a_over_2_v3_deriv(1:3, 1:3, nb) = director_deriv(1:3, 1:3)
785  a_over_2_v3_second(1:3, 1:3, 1:3, nb) = director_second(1:3, 1:3, 1:3)
786  endif
787  endif
788  if( .not. use_director_tangent ) then
789  a_over_2_v3_deriv(1:3, 1, nb) = (/ 0.0d0, -a_over_2_v3(3, nb), a_over_2_v3(2, nb) /)
790  a_over_2_v3_deriv(1:3, 2, nb) = (/ a_over_2_v3(3, nb), 0.0d0, -a_over_2_v3(1, nb) /)
791  a_over_2_v3_deriv(1:3, 3, nb) = (/ -a_over_2_v3(2, nb), a_over_2_v3(1, nb), 0.0d0 /)
792  endif
793 
794  !--------------------------------------------------------
795 
796  end do
797 
798  !--------------------------------------------------------------------
799  ! MODIFIED to LAMINATED SHELL ANALYSIS
800  !--------------------------------------------------------------------
801  zeta_ly = 0.0d0
802 
803  n_totlyr = gausses(1)%pMaterial%totallyr
804  do n_layer=1,n_totlyr
805  do ly = 1, ny
806 
807  !--------------------------------------------------------
808 
809  ! MITC4
810  if( etype .EQ. fe_mitc4_shell ) then
811 
812  zeta_ly = 0.0d0
813 
814  ! MITC9
815  else if( etype .EQ. fe_mitc9_shell ) then
816 
817  zeta_ly = gauss1d3(1,ly)
818 
819  ! MITC3
820  else if( etype .EQ. fe_mitc3_shell )then
821 
822  zeta_ly = 0.0d0
823 
824  end if
825 
826  !--------------------------------------------------------
827 
828  do it = 1, ntying
829 
830  do ip = 1, npoints_tying(it)
831 
832  !-------------------------------------------------
833 
834  naturalcoord(1) = tpcoord(ip, 1, it)
835  naturalcoord(2) = tpcoord(ip, 2, it)
836 
837  call getshapefunc(fetype, naturalcoord, shapefunc)
838 
839  call getshapederiv(fetype, naturalcoord, shapederiv)
840 
841  !-------------------------------------------------
842 
843  do na = 1, nn
844 
845  do i = 1, 3
846 
847  u_rot(i, na) &
848  = shapefunc(na) &
849  *( zeta_ly*a_over_2_v3(i, na) )
850 
851  dudxi_rot(i, na) &
852  = shapederiv(na, 1) &
853  *( zeta_ly*a_over_2_v3(i, na) )
854  dudeta_rot(i, na) &
855  = shapederiv(na, 2) &
856  *( zeta_ly*a_over_2_v3(i, na) )
857  dudzeta_rot(i, na) &
858  = shapefunc(na) &
859  *( a_over_2_v3(i, na) )
860 
861  end do
862 
863  end do
864 
865  !-------------------------------------------------
866 
867  ! Covariant basis vector
868  do i = 1, 3
869 
870  g1(i) = 0.0d0
871  g2(i) = 0.0d0
872  g3(i) = 0.0d0
873 
874  do na = 1, nn
875 
876  g1(i) = g1(i)+shapederiv(na, 1) &
877  *elem(i, na) &
878  +dudxi_rot(i, na)
879  g2(i) = g2(i)+shapederiv(na, 2) &
880  *elem(i, na) &
881  +dudeta_rot(i, na)
882  g3(i) = g3(i)+dudzeta_rot(i, na)
883 
884  end do
885 
886  end do
887 
888  !-------------------------------------------------
889 
890  g1_tl(:) = g1(:)
891  g2_tl(:) = g2(:)
892  if( use_tl_green ) then
893  dudxi_trans(:) = 0.0d0
894  dudeta_trans(:) = 0.0d0
895  do na = 1, nn
896  dudxi_trans(:) = dudxi_trans(:)+shapederiv(na, 1)*shell_disp(1:3, na)
897  dudeta_trans(:) = dudeta_trans(:)+shapederiv(na, 2)*shell_disp(1:3, na)
898  end do
899  if( etype == fe_mitc9_shell ) then
900  g1_tl(:) = g1_tl(:)+dudxi_trans(:)
901  g2_tl(:) = g2_tl(:)+dudeta_trans(:)
902  else
903  g1(:) = g1(:)+dudxi_trans(:)
904  g2(:) = g2(:)+dudeta_trans(:)
905  g1_tl(:) = g1(:)
906  g2_tl(:) = g2(:)
907  endif
908  endif
909 
910  !-------------------------------------------------
911 
912  ! [ B L ] matrix
913  do nb = 1, nn
914 
915  jsize1 = ndof*(nb-1)+1
916  jsize2 = ndof*(nb-1)+2
917  jsize3 = ndof*(nb-1)+3
918  jsize4 = ndof*(nb-1)+4
919  jsize5 = ndof*(nb-1)+5
920  jsize6 = ndof*(nb-1)+6
921 
922  aa1(1) = dudxi_rot(2, nb) *g1(3)-dudxi_rot(3, nb) *g1(2)
923  aa1(2) = dudxi_rot(3, nb) *g1(1)-dudxi_rot(1, nb) *g1(3)
924  aa1(3) = dudxi_rot(1, nb) *g1(2)-dudxi_rot(2, nb) *g1(1)
925 
926  aa2(1) = dudxi_rot(2, nb) *g2(3)-dudxi_rot(3, nb) *g2(2)
927  aa2(2) = dudxi_rot(3, nb) *g2(1)-dudxi_rot(1, nb) *g2(3)
928  aa2(3) = dudxi_rot(1, nb) *g2(2)-dudxi_rot(2, nb) *g2(1)
929 
930  aa3(1) = dudxi_rot(2, nb) *g3(3)-dudxi_rot(3, nb) *g3(2)
931  aa3(2) = dudxi_rot(3, nb) *g3(1)-dudxi_rot(1, nb) *g3(3)
932  aa3(3) = dudxi_rot(1, nb) *g3(2)-dudxi_rot(2, nb) *g3(1)
933 
934  bb1(1) = dudeta_rot(2, nb) *g1(3)-dudeta_rot(3, nb) *g1(2)
935  bb1(2) = dudeta_rot(3, nb) *g1(1)-dudeta_rot(1, nb) *g1(3)
936  bb1(3) = dudeta_rot(1, nb) *g1(2)-dudeta_rot(2, nb) *g1(1)
937 
938  bb2(1) = dudeta_rot(2, nb) *g2(3)-dudeta_rot(3, nb) *g2(2)
939  bb2(2) = dudeta_rot(3, nb) *g2(1)-dudeta_rot(1, nb) *g2(3)
940  bb2(3) = dudeta_rot(1, nb) *g2(2)-dudeta_rot(2, nb) *g2(1)
941 
942  bb3(1) = dudeta_rot(2, nb) *g3(3)-dudeta_rot(3, nb) *g3(2)
943  bb3(2) = dudeta_rot(3, nb) *g3(1)-dudeta_rot(1, nb) *g3(3)
944  bb3(3) = dudeta_rot(1, nb) *g3(2)-dudeta_rot(2, nb) *g3(1)
945 
946  cc1(1) = dudzeta_rot(2, nb)*g1(3)-dudzeta_rot(3, nb)*g1(2)
947  cc1(2) = dudzeta_rot(3, nb)*g1(1)-dudzeta_rot(1, nb)*g1(3)
948  cc1(3) = dudzeta_rot(1, nb)*g1(2)-dudzeta_rot(2, nb)*g1(1)
949 
950  cc2(1) = dudzeta_rot(2, nb)*g2(3)-dudzeta_rot(3, nb)*g2(2)
951  cc2(2) = dudzeta_rot(3, nb)*g2(1)-dudzeta_rot(1, nb)*g2(3)
952  cc2(3) = dudzeta_rot(1, nb)*g2(2)-dudzeta_rot(2, nb)*g2(1)
953 
954  if( use_director_tangent ) then
955  dudxi_rot_deriv(1:3, 1:3) = shapederiv(nb, 1)*zeta_ly &
956  *a_over_2_v3_deriv(1:3, 1:3, nb)
957  dudeta_rot_deriv(1:3, 1:3) = shapederiv(nb, 2)*zeta_ly &
958  *a_over_2_v3_deriv(1:3, 1:3, nb)
959  dudzeta_rot_deriv(1:3, 1:3) = shapefunc(nb) &
960  *a_over_2_v3_deriv(1:3, 1:3, nb)
961  dudxi_rot_second(1:3, 1:3, 1:3) = shapederiv(nb, 1)*zeta_ly &
962  *a_over_2_v3_second(1:3, 1:3, 1:3, nb)
963  dudeta_rot_second(1:3, 1:3, 1:3) = shapederiv(nb, 2)*zeta_ly &
964  *a_over_2_v3_second(1:3, 1:3, 1:3, nb)
965  dudzeta_rot_second(1:3, 1:3, 1:3) = shapefunc(nb) &
966  *a_over_2_v3_second(1:3, 1:3, 1:3, nb)
967  endif
968 
969  b_di(1, jsize1, ip, it, ly) = shapederiv(nb, 1)*g1_tl(1)
970  b_di(2, jsize1, ip, it, ly) = shapederiv(nb, 2)*g2_tl(1)
971  b_di(3, jsize1, ip, it, ly) = shapederiv(nb, 1)*g2_tl(1) &
972  +shapederiv(nb, 2)*g1_tl(1)
973  b_di(4, jsize1, ip, it, ly) = shapederiv(nb, 2)*g3(1)
974  b_di(5, jsize1, ip, it, ly) = shapederiv(nb, 1)*g3(1)
975 
976  b_di(1, jsize2, ip, it, ly) = shapederiv(nb, 1)*g1_tl(2)
977  b_di(2, jsize2, ip, it, ly) = shapederiv(nb, 2)*g2_tl(2)
978  b_di(3, jsize2, ip, it, ly) = shapederiv(nb, 1)*g2_tl(2) &
979  +shapederiv(nb, 2)*g1_tl(2)
980  b_di(4, jsize2, ip, it, ly) = shapederiv(nb, 2)*g3(2)
981  b_di(5, jsize2, ip, it, ly) = shapederiv(nb, 1)*g3(2)
982 
983  b_di(1, jsize3, ip, it, ly) = shapederiv(nb, 1)*g1_tl(3)
984  b_di(2, jsize3, ip, it, ly) = shapederiv(nb, 2)*g2_tl(3)
985  b_di(3, jsize3, ip, it, ly) = shapederiv(nb, 1)*g2_tl(3) &
986  +shapederiv(nb, 2)*g1_tl(3)
987  b_di(4, jsize3, ip, it, ly) = shapederiv(nb, 2)*g3(3)
988  b_di(5, jsize3, ip, it, ly) = shapederiv(nb, 1)*g3(3)
989 
990  if( use_director_tangent ) then
991  do m = 1, 3
992  jsize = ndof*(nb-1)+3+m
993  b_di(1, jsize, ip, it, ly) = dot_product(dudxi_rot_deriv(1:3, m), g1)
994  b_di(2, jsize, ip, it, ly) = dot_product(dudeta_rot_deriv(1:3, m), g2)
995  b_di(3, jsize, ip, it, ly) = dot_product(dudxi_rot_deriv(1:3, m), g2) &
996  +dot_product(dudeta_rot_deriv(1:3, m), g1)
997  b_di(4, jsize, ip, it, ly) = dot_product(dudeta_rot_deriv(1:3, m), g3) &
998  +dot_product(dudzeta_rot_deriv(1:3, m), g2)
999  b_di(5, jsize, ip, it, ly) = dot_product(dudxi_rot_deriv(1:3, m), g3) &
1000  +dot_product(dudzeta_rot_deriv(1:3, m), g1)
1001  end do
1002  do n = 1, 3
1003  do m = 1, 3
1004  b2rot_di(1, m, n, nb, ip, it, ly) = dot_product(dudxi_rot_second(1:3, m, n), g1)
1005  b2rot_di(2, m, n, nb, ip, it, ly) = dot_product(dudeta_rot_second(1:3, m, n), g2)
1006  b2rot_di(3, m, n, nb, ip, it, ly) = dot_product(dudxi_rot_second(1:3, m, n), g2) &
1007  +dot_product(dudeta_rot_second(1:3, m, n), g1)
1008  b2rot_di(4, m, n, nb, ip, it, ly) = dot_product(dudeta_rot_second(1:3, m, n), g3) &
1009  +dot_product(dudzeta_rot_second(1:3, m, n), g2)
1010  b2rot_di(5, m, n, nb, ip, it, ly) = dot_product(dudxi_rot_second(1:3, m, n), g3) &
1011  +dot_product(dudzeta_rot_second(1:3, m, n), g1)
1012  end do
1013  end do
1014  else
1015  b_di(1, jsize4, ip, it, ly) = aa1(1)
1016  b_di(2, jsize4, ip, it, ly) = bb2(1)
1017  b_di(3, jsize4, ip, it, ly) = aa2(1)+bb1(1)
1018  b_di(4, jsize4, ip, it, ly) = bb3(1)+cc2(1)
1019  b_di(5, jsize4, ip, it, ly) = aa3(1)+cc1(1)
1020 
1021  b_di(1, jsize5, ip, it, ly) = aa1(2)
1022  b_di(2, jsize5, ip, it, ly) = bb2(2)
1023  b_di(3, jsize5, ip, it, ly) = aa2(2)+bb1(2)
1024  b_di(4, jsize5, ip, it, ly) = bb3(2)+cc2(2)
1025  b_di(5, jsize5, ip, it, ly) = aa3(2)+cc1(2)
1026 
1027  b_di(1, jsize6, ip, it, ly) = aa1(3)
1028  b_di(2, jsize6, ip, it, ly) = bb2(3)
1029  b_di(3, jsize6, ip, it, ly) = aa2(3)+bb1(3)
1030  b_di(4, jsize6, ip, it, ly) = bb3(3)+cc2(3)
1031  b_di(5, jsize6, ip, it, ly) = aa3(3)+cc1(3)
1032  endif
1033 
1034  ! First variations of g1/g2/g3 at tying points, used by the MITC4 shear geometric stiffness.
1035  bg1_di(1:3, jsize1, ip, it, ly) = (/ shapederiv(nb, 1), 0.0d0, 0.0d0 /)
1036  bg1_di(1:3, jsize2, ip, it, ly) = (/ 0.0d0, shapederiv(nb, 1), 0.0d0 /)
1037  bg1_di(1:3, jsize3, ip, it, ly) = (/ 0.0d0, 0.0d0, shapederiv(nb, 1) /)
1038  bg2_di(1:3, jsize1, ip, it, ly) = (/ shapederiv(nb, 2), 0.0d0, 0.0d0 /)
1039  bg2_di(1:3, jsize2, ip, it, ly) = (/ 0.0d0, shapederiv(nb, 2), 0.0d0 /)
1040  bg2_di(1:3, jsize3, ip, it, ly) = (/ 0.0d0, 0.0d0, shapederiv(nb, 2) /)
1041  if( use_director_tangent ) then
1042  do m = 1, 3
1043  jsize = ndof*(nb-1)+3+m
1044  bg1_di(1:3, jsize, ip, it, ly) = dudxi_rot_deriv(1:3, m)
1045  bg2_di(1:3, jsize, ip, it, ly) = dudeta_rot_deriv(1:3, m)
1046  bg3_di(1:3, jsize, ip, it, ly) = dudzeta_rot_deriv(1:3, m)
1047  end do
1048  else
1049  bg1_di(1:3, jsize4, ip, it, ly) = (/ 0.0d0, -dudxi_rot(3, nb), dudxi_rot(2, nb) /)
1050  bg1_di(1:3, jsize5, ip, it, ly) = (/ dudxi_rot(3, nb), 0.0d0, -dudxi_rot(1, nb) /)
1051  bg1_di(1:3, jsize6, ip, it, ly) = (/ -dudxi_rot(2, nb), dudxi_rot(1, nb), 0.0d0 /)
1052  bg2_di(1:3, jsize4, ip, it, ly) = (/ 0.0d0, -dudeta_rot(3, nb), dudeta_rot(2, nb) /)
1053  bg2_di(1:3, jsize5, ip, it, ly) = (/ dudeta_rot(3, nb), 0.0d0, -dudeta_rot(1, nb) /)
1054  bg2_di(1:3, jsize6, ip, it, ly) = (/ -dudeta_rot(2, nb), dudeta_rot(1, nb), 0.0d0 /)
1055  bg3_di(1:3, jsize4, ip, it, ly) = (/ 0.0d0, -dudzeta_rot(3, nb), dudzeta_rot(2, nb) /)
1056  bg3_di(1:3, jsize5, ip, it, ly) = (/ dudzeta_rot(3, nb), 0.0d0, -dudzeta_rot(1, nb) /)
1057  bg3_di(1:3, jsize6, ip, it, ly) = (/ -dudzeta_rot(2, nb), dudzeta_rot(1, nb), 0.0d0 /)
1058  endif
1059 
1060  end do
1061 
1062  !-------------------------------------------------
1063 
1064  end do
1065 
1066  end do
1067 
1068  !--------------------------------------------------------
1069 
1070  call fstr_shell_layer_quadrature_gauss( etype, gausses(1), n_layer, ly, &
1071  zeta_ly, w_ly, ierr_quad )
1072  if( ierr_quad /= 0 ) cycle
1073 
1074  !--------------------------------------------------------
1075 
1076  do lx = 1, numofquadpoints(fetype)
1077 
1078  !--------------------------------------------------
1079 
1080  call getquadpoint(fetype, lx, naturalcoord)
1081 
1082  xi_lx = naturalcoord(1)
1083  eta_lx = naturalcoord(2)
1084 
1085  w_w_lx = getweight(fetype, lx)
1086 
1087  call getshapefunc(fetype, naturalcoord, shapefunc)
1088 
1089  call getshapederiv(fetype, naturalcoord, shapederiv)
1090 
1091  !--------------------------------------------------
1092 
1093  do i = 1, 3
1094 
1095  v1_i(i) = 0.0d0
1096  v2_i(i) = 0.0d0
1097  v3_i(i) = 0.0d0
1098 
1099  do na = 1, nn
1100 
1101  v1_i(i) = v1_i(i)+shapefunc(na)*v1(i, na)
1102  v2_i(i) = v2_i(i)+shapefunc(na)*v2(i, na)
1103  v3_i(i) = v3_i(i)+shapefunc(na)*v3(i, na)
1104 
1105  end do
1106 
1107  end do
1108 
1109  !--------------------------------------------------
1110 
1111  do na = 1, nn
1112 
1113  do i = 1, 3
1114 
1115  u_rot(i, na) &
1116  = shapefunc(na) &
1117  *( zeta_ly*a_over_2_v3(i, na) )
1118 
1119  dudxi_rot(i, na) &
1120  = shapederiv(na, 1) &
1121  *( zeta_ly*a_over_2_v3(i, na) )
1122  dudeta_rot(i, na) &
1123  = shapederiv(na, 2) &
1124  *( zeta_ly*a_over_2_v3(i, na) )
1125  dudzeta_rot(i, na) &
1126  = shapefunc(na) &
1127  *( a_over_2_v3(i, na) )
1128 
1129  end do
1130 
1131  end do
1132 
1133  !--------------------------------------------------
1134 
1135  ! Covariant basis vector
1136  do i = 1, 3
1137 
1138  g1(i) = 0.0d0
1139  g2(i) = 0.0d0
1140  g3(i) = 0.0d0
1141 
1142  do na = 1, nn
1143 
1144  g1(i) = g1(i)+shapederiv(na, 1) &
1145  *elem(i, na) &
1146  +dudxi_rot(i, na)
1147  g2(i) = g2(i)+shapederiv(na, 2) &
1148  *elem(i, na) &
1149  +dudeta_rot(i, na)
1150  g3(i) = g3(i)+dudzeta_rot(i, na)
1151 
1152  end do
1153 
1154  end do
1155 
1156  !--------------------------------------------------
1157 
1158  if( use_tl_green ) then
1159  ! reference shell geometry for TL quadrature weight
1160  g1_weight(:) = 0.0d0
1161  g2_weight(:) = 0.0d0
1162  g3_weight(:) = 0.0d0
1163  do na = 1, nn
1164  g1_weight(:) = g1_weight(:)+shapederiv(na, 1) &
1165  *(ecoord(1:3, na)+zeta_ly*a_over_2_v3_ref(1:3, na))
1166  g2_weight(:) = g2_weight(:)+shapederiv(na, 2) &
1167  *(ecoord(1:3, na)+zeta_ly*a_over_2_v3_ref(1:3, na))
1168  g3_weight(:) = g3_weight(:)+shapefunc(na)*a_over_2_v3_ref(1:3, na)
1169  end do
1170  det_weight = g1_weight(1)*( g2_weight(2)*g3_weight(3)-g2_weight(3)*g3_weight(2) ) &
1171  +g1_weight(2)*( g2_weight(3)*g3_weight(1)-g2_weight(1)*g3_weight(3) ) &
1172  +g1_weight(3)*( g2_weight(1)*g3_weight(2)-g2_weight(2)*g3_weight(1) )
1173  endif
1174 
1175  !--------------------------------------------------
1176 
1177  if( use_tl_green ) then
1178  ! Green-Lagrange membrane and bending strain
1179  dudxi_trans(:) = 0.0d0
1180  dudeta_trans(:) = 0.0d0
1181  do na = 1, nn
1182  dudxi_trans(:) = dudxi_trans(:)+shapederiv(na, 1)*shell_disp(1:3, na)
1183  dudeta_trans(:) = dudeta_trans(:)+shapederiv(na, 2)*shell_disp(1:3, na)
1184  end do
1185  g1(:) = g1(:)+dudxi_trans(:)
1186  g2(:) = g2(:)+dudeta_trans(:)
1187  endif
1188 
1189  !--------------------------------------------------
1190 
1191  ! Jacobian
1192  det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
1193  +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
1194  +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )
1195 
1196  if( .not. use_tl_green ) det_weight = det
1197 
1198  det_inv = 1.0d0/det
1199 
1200  !--------------------------------------------------
1201 
1202  ! Contravariant basis vector
1203  cg1(1) = det_inv &
1204  *( g2(2)*g3(3)-g2(3)*g3(2) )
1205  cg1(2) = det_inv &
1206  *( g2(3)*g3(1)-g2(1)*g3(3) )
1207  cg1(3) = det_inv &
1208  *( g2(1)*g3(2)-g2(2)*g3(1) )
1209  cg2(1) = det_inv &
1210  *( g3(2)*g1(3)-g3(3)*g1(2) )
1211  cg2(2) = det_inv &
1212  *( g3(3)*g1(1)-g3(1)*g1(3) )
1213  cg2(3) = det_inv &
1214  *( g3(1)*g1(2)-g3(2)*g1(1) )
1215  cg3(1) = det_inv &
1216  *( g1(2)*g2(3)-g1(3)*g2(2) )
1217  cg3(2) = det_inv &
1218  *( g1(3)*g2(1)-g1(1)*g2(3) )
1219  cg3(3) = det_inv &
1220  *( g1(1)*g2(2)-g1(2)*g2(1) )
1221 
1222  !--------------------------------------------------
1223 
1224  g3_abs = dsqrt( g3(1)*g3(1) &
1225  +g3(2)*g3(2) &
1226  +g3(3)*g3(3) )
1227 
1228  !--------------------------------------------------
1229 
1230  ! Orthonormal vectors
1231 
1232  e3_hat(1) = g3(1)/g3_abs
1233  e3_hat(2) = g3(2)/g3_abs
1234  e3_hat(3) = g3(3)/g3_abs
1235 
1236  e1_hat(1) = g2(2)*e3_hat(3) &
1237  -g2(3)*e3_hat(2)
1238  e1_hat(2) = g2(3)*e3_hat(1) &
1239  -g2(1)*e3_hat(3)
1240  e1_hat(3) = g2(1)*e3_hat(2) &
1241  -g2(2)*e3_hat(1)
1242  e1_hat_abs = dsqrt( e1_hat(1)*e1_hat(1) &
1243  +e1_hat(2)*e1_hat(2) &
1244  +e1_hat(3)*e1_hat(3) )
1245  e1_hat(1) = e1_hat(1)/e1_hat_abs
1246  e1_hat(2) = e1_hat(2)/e1_hat_abs
1247  e1_hat(3) = e1_hat(3)/e1_hat_abs
1248 
1249  e2_hat(1) = e3_hat(2)*e1_hat(3) &
1250  -e3_hat(3)*e1_hat(2)
1251  e2_hat(2) = e3_hat(3)*e1_hat(1) &
1252  -e3_hat(1)*e1_hat(3)
1253  e2_hat(3) = e3_hat(1)*e1_hat(2) &
1254  -e3_hat(2)*e1_hat(1)
1255  e2_hat_abs = dsqrt( e2_hat(1)*e2_hat(1) &
1256  +e2_hat(2)*e2_hat(2) &
1257  +e2_hat(3)*e2_hat(3) )
1258  e2_hat(1) = e2_hat(1)/e2_hat_abs
1259  e2_hat(2) = e2_hat(2)/e2_hat_abs
1260  e2_hat(3) = e2_hat(3)/e2_hat_abs
1261 
1262  !--------------------------------------------------
1263 
1264  e1_hat_mat(:) = e1_hat(:)
1265  e2_hat_mat(:) = e2_hat(:)
1266  e3_hat_mat(:) = e3_hat(:)
1267  cg1_mat(:) = cg1(:)
1268  cg2_mat(:) = cg2(:)
1269  cg3_mat(:) = cg3(:)
1270  if( use_tl_green ) then
1271  call shell_basis_from_covariant(g1_weight, g2_weight, g3_weight, &
1272  e1_hat_mat, e2_hat_mat, e3_hat_mat, cg1_mat, cg2_mat, cg3_mat, det_mat)
1273  endif
1274 
1275  !--------------------------------------------------
1276 
1277  if( present( element ) ) then
1278  ishell = fstr_shell_layer_gauss_index( element, lx, n_layer, ly )
1279  else
1280  ishell = 0
1281  endif
1282 
1283  if( ishell > 0 ) then
1284  call matlmatrix_shell &
1285  (element%shell_layer_gausses(ishell), shell, d, &
1286  e1_hat_mat, e2_hat_mat, e3_hat_mat, cg1_mat, cg2_mat, cg3_mat, &
1287  alpha, n_layer)
1288  else
1289  call matlmatrix_shell &
1290  (gausses(lx), shell, d, &
1291  e1_hat_mat, e2_hat_mat, e3_hat_mat, cg1_mat, cg2_mat, cg3_mat, &
1292  alpha, n_layer)
1293  endif
1294  !--------------------------------------------------
1295 
1296  g1_tl(:) = g1(:)
1297  g2_tl(:) = g2(:)
1298  if( use_tl_green ) then
1299  dudxi_trans(:) = 0.0d0
1300  dudeta_trans(:) = 0.0d0
1301  do na = 1, nn
1302  dudxi_trans(:) = dudxi_trans(:)+shapederiv(na, 1)*shell_disp(1:3, na)
1303  dudeta_trans(:) = dudeta_trans(:)+shapederiv(na, 2)*shell_disp(1:3, na)
1304  end do
1305  if( etype == fe_mitc9_shell ) then
1306  g1_tl(:) = g1_tl(:)+dudxi_trans(:)
1307  g2_tl(:) = g2_tl(:)+dudeta_trans(:)
1308  else if( etype /= fe_mitc4_shell ) then
1309  g1(:) = g1(:)+dudxi_trans(:)
1310  g2(:) = g2(:)+dudeta_trans(:)
1311  g1_tl(:) = g1(:)
1312  g2_tl(:) = g2(:)
1313  endif
1314  endif
1315 
1316  !--------------------------------------------------
1317 
1318  ! [ B L ] matrix
1319  b(:, :) = 0.0d0
1320  if( use_director_tangent ) b2rot(:, :, :, :) = 0.0d0
1321  do nb = 1, nn
1322 
1323  jsize1 = ndof*(nb-1)+1
1324  jsize2 = ndof*(nb-1)+2
1325  jsize3 = ndof*(nb-1)+3
1326  jsize4 = ndof*(nb-1)+4
1327  jsize5 = ndof*(nb-1)+5
1328  jsize6 = ndof*(nb-1)+6
1329 
1330  aa1(1) = dudxi_rot(2, nb) *g1(3)-dudxi_rot(3, nb) *g1(2)
1331  aa1(2) = dudxi_rot(3, nb) *g1(1)-dudxi_rot(1, nb) *g1(3)
1332  aa1(3) = dudxi_rot(1, nb) *g1(2)-dudxi_rot(2, nb) *g1(1)
1333 
1334  aa2(1) = dudxi_rot(2, nb) *g2(3)-dudxi_rot(3, nb) *g2(2)
1335  aa2(2) = dudxi_rot(3, nb) *g2(1)-dudxi_rot(1, nb) *g2(3)
1336  aa2(3) = dudxi_rot(1, nb) *g2(2)-dudxi_rot(2, nb) *g2(1)
1337 
1338  aa3(1) = dudxi_rot(2, nb) *g3(3)-dudxi_rot(3, nb) *g3(2)
1339  aa3(2) = dudxi_rot(3, nb) *g3(1)-dudxi_rot(1, nb) *g3(3)
1340  aa3(3) = dudxi_rot(1, nb) *g3(2)-dudxi_rot(2, nb) *g3(1)
1341 
1342  bb1(1) = dudeta_rot(2, nb) *g1(3)-dudeta_rot(3, nb) *g1(2)
1343  bb1(2) = dudeta_rot(3, nb) *g1(1)-dudeta_rot(1, nb) *g1(3)
1344  bb1(3) = dudeta_rot(1, nb) *g1(2)-dudeta_rot(2, nb) *g1(1)
1345 
1346  bb2(1) = dudeta_rot(2, nb) *g2(3)-dudeta_rot(3, nb) *g2(2)
1347  bb2(2) = dudeta_rot(3, nb) *g2(1)-dudeta_rot(1, nb) *g2(3)
1348  bb2(3) = dudeta_rot(1, nb) *g2(2)-dudeta_rot(2, nb) *g2(1)
1349 
1350  bb3(1) = dudeta_rot(2, nb) *g3(3)-dudeta_rot(3, nb) *g3(2)
1351  bb3(2) = dudeta_rot(3, nb) *g3(1)-dudeta_rot(1, nb) *g3(3)
1352  bb3(3) = dudeta_rot(1, nb) *g3(2)-dudeta_rot(2, nb) *g3(1)
1353 
1354  cc1(1) = dudzeta_rot(2, nb)*g1(3)-dudzeta_rot(3, nb)*g1(2)
1355  cc1(2) = dudzeta_rot(3, nb)*g1(1)-dudzeta_rot(1, nb)*g1(3)
1356  cc1(3) = dudzeta_rot(1, nb)*g1(2)-dudzeta_rot(2, nb)*g1(1)
1357 
1358  cc2(1) = dudzeta_rot(2, nb)*g2(3)-dudzeta_rot(3, nb)*g2(2)
1359  cc2(2) = dudzeta_rot(3, nb)*g2(1)-dudzeta_rot(1, nb)*g2(3)
1360  cc2(3) = dudzeta_rot(1, nb)*g2(2)-dudzeta_rot(2, nb)*g2(1)
1361 
1362  if( use_director_tangent ) then
1363  dudxi_rot_deriv(1:3, 1:3) = shapederiv(nb, 1)*zeta_ly &
1364  *a_over_2_v3_deriv(1:3, 1:3, nb)
1365  dudeta_rot_deriv(1:3, 1:3) = shapederiv(nb, 2)*zeta_ly &
1366  *a_over_2_v3_deriv(1:3, 1:3, nb)
1367  dudzeta_rot_deriv(1:3, 1:3) = shapefunc(nb) &
1368  *a_over_2_v3_deriv(1:3, 1:3, nb)
1369  dudxi_rot_second(1:3, 1:3, 1:3) = shapederiv(nb, 1)*zeta_ly &
1370  *a_over_2_v3_second(1:3, 1:3, 1:3, nb)
1371  dudeta_rot_second(1:3, 1:3, 1:3) = shapederiv(nb, 2)*zeta_ly &
1372  *a_over_2_v3_second(1:3, 1:3, 1:3, nb)
1373  dudzeta_rot_second(1:3, 1:3, 1:3) = shapefunc(nb) &
1374  *a_over_2_v3_second(1:3, 1:3, 1:3, nb)
1375  endif
1376 
1377  b(1, jsize1) = shapederiv(nb, 1)*g1_tl(1)
1378  b(2, jsize1) = shapederiv(nb, 2)*g2_tl(1)
1379  b(3, jsize1) = shapederiv(nb, 1)*g2_tl(1) &
1380  +shapederiv(nb, 2)*g1_tl(1)
1381  b(4, jsize1) = shapederiv(nb, 2)*g3(1)
1382  b(5, jsize1) = shapederiv(nb, 1)*g3(1)
1383 
1384  b(1, jsize2) = shapederiv(nb, 1)*g1_tl(2)
1385  b(2, jsize2) = shapederiv(nb, 2)*g2_tl(2)
1386  b(3, jsize2) = shapederiv(nb, 1)*g2_tl(2) &
1387  +shapederiv(nb, 2)*g1_tl(2)
1388  b(4, jsize2) = shapederiv(nb, 2)*g3(2)
1389  b(5, jsize2) = shapederiv(nb, 1)*g3(2)
1390 
1391  b(1, jsize3) = shapederiv(nb, 1)*g1_tl(3)
1392  b(2, jsize3) = shapederiv(nb, 2)*g2_tl(3)
1393  b(3, jsize3) = shapederiv(nb, 1)*g2_tl(3) &
1394  +shapederiv(nb, 2)*g1_tl(3)
1395  b(4, jsize3) = shapederiv(nb, 2)*g3(3)
1396  b(5, jsize3) = shapederiv(nb, 1)*g3(3)
1397 
1398  if( use_director_tangent ) then
1399  do m = 1, 3
1400  jsize = ndof*(nb-1)+3+m
1401  b(1, jsize) = dot_product(dudxi_rot_deriv(1:3, m), g1)
1402  b(2, jsize) = dot_product(dudeta_rot_deriv(1:3, m), g2)
1403  b(3, jsize) = dot_product(dudxi_rot_deriv(1:3, m), g2) &
1404  +dot_product(dudeta_rot_deriv(1:3, m), g1)
1405  b(4, jsize) = dot_product(dudeta_rot_deriv(1:3, m), g3) &
1406  +dot_product(dudzeta_rot_deriv(1:3, m), g2)
1407  b(5, jsize) = dot_product(dudxi_rot_deriv(1:3, m), g3) &
1408  +dot_product(dudzeta_rot_deriv(1:3, m), g1)
1409  end do
1410  do n = 1, 3
1411  do m = 1, 3
1412  b2rot(1, m, n, nb) = dot_product(dudxi_rot_second(1:3, m, n), g1)
1413  b2rot(2, m, n, nb) = dot_product(dudeta_rot_second(1:3, m, n), g2)
1414  b2rot(3, m, n, nb) = dot_product(dudxi_rot_second(1:3, m, n), g2) &
1415  +dot_product(dudeta_rot_second(1:3, m, n), g1)
1416  b2rot(4, m, n, nb) = dot_product(dudeta_rot_second(1:3, m, n), g3) &
1417  +dot_product(dudzeta_rot_second(1:3, m, n), g2)
1418  b2rot(5, m, n, nb) = dot_product(dudxi_rot_second(1:3, m, n), g3) &
1419  +dot_product(dudzeta_rot_second(1:3, m, n), g1)
1420  end do
1421  end do
1422  else
1423  b(1, jsize4) = aa1(1)
1424  b(2, jsize4) = bb2(1)
1425  b(3, jsize4) = aa2(1)+bb1(1)
1426  b(4, jsize4) = bb3(1)+cc2(1)
1427  b(5, jsize4) = aa3(1)+cc1(1)
1428 
1429  b(1, jsize5) = aa1(2)
1430  b(2, jsize5) = bb2(2)
1431  b(3, jsize5) = aa2(2)+bb1(2)
1432  b(4, jsize5) = bb3(2)+cc2(2)
1433  b(5, jsize5) = aa3(2)+cc1(2)
1434 
1435  b(1, jsize6) = aa1(3)
1436  b(2, jsize6) = bb2(3)
1437  b(3, jsize6) = aa2(3)+bb1(3)
1438  b(4, jsize6) = bb3(3)+cc2(3)
1439  b(5, jsize6) = aa3(3)+cc1(3)
1440  endif
1441 
1442  end do
1443 
1444  !--------------------------------------------------
1445 
1446  ! MITC4
1447  if( etype .EQ. fe_mitc4_shell ) then
1448 
1449  do jsize = 1, ndof*nn
1450 
1451  b(4, jsize) = 0.0d0
1452  b(5, jsize) = 0.0d0
1453 
1454  ! B_as(4, jsize)
1455  b(4, jsize) &
1456  = 0.5d0*( 1.0d0-xi_lx )*b_di(4, jsize, 4, 1, ly) &
1457  +0.5d0*( 1.0d0+xi_lx )*b_di(4, jsize, 2, 1, ly)
1458  ! B_as(5, jsize)
1459  b(5, jsize) &
1460  = 0.5d0*( 1.0d0-eta_lx )*b_di(5, jsize, 1, 1, ly) &
1461  +0.5d0*( 1.0d0+eta_lx )*b_di(5, jsize, 3, 1, ly)
1462 
1463  end do
1464  if( use_director_tangent ) then
1465  do nb = 1, nn
1466  do n = 1, 3
1467  do m = 1, 3
1468  b2rot(4, m, n, nb) &
1469  = 0.5d0*( 1.0d0-xi_lx )*b2rot_di(4, m, n, nb, 4, 1, ly) &
1470  +0.5d0*( 1.0d0+xi_lx )*b2rot_di(4, m, n, nb, 2, 1, ly)
1471  b2rot(5, m, n, nb) &
1472  = 0.5d0*( 1.0d0-eta_lx )*b2rot_di(5, m, n, nb, 1, 1, ly) &
1473  +0.5d0*( 1.0d0+eta_lx )*b2rot_di(5, m, n, nb, 3, 1, ly)
1474  end do
1475  end do
1476  end do
1477  endif
1478 
1479  ! MITC9
1480  else if( etype .EQ. fe_mitc9_shell ) then
1481 
1482  xxi_lx = xi_lx /dsqrt( 1.0d0/3.0d0 )
1483  eeta_lx = eta_lx/dsqrt( 3.0d0/5.0d0 )
1484 
1485  do ip = 1, npoints_tying(1)
1486 
1487  h(ip, 1) &
1488  = ( 0.5d0*( 1.0d0+xxi_di(ip, 1)*xxi_lx ) ) &
1489  *( ( 0.5d0*eeta_di(ip, 1)*eeta_lx ) &
1490  *( 1.0d0+eeta_di(ip, 1)*eeta_lx ) &
1491  +( 1.0d0-eeta_di(ip, 1)*eeta_di(ip, 1) ) &
1492  *( 1.0d0-eeta_lx*eeta_lx ) )
1493 
1494  end do
1495 
1496  xxi_lx = xi_lx /dsqrt( 3.0d0/5.0d0 )
1497  eeta_lx = eta_lx/dsqrt( 1.0d0/3.0d0 )
1498 
1499  do ip = 1, npoints_tying(2)
1500 
1501  h(ip, 2) &
1502  = ( ( 0.5d0*xxi_di(ip, 2) *xxi_lx ) &
1503  *( 1.0d0+xxi_di(ip, 2) *xxi_lx ) &
1504  +( 1.0d0-xxi_di(ip, 2) *xxi_di(ip, 2) ) &
1505  *( 1.0d0-xxi_lx*xxi_lx ) ) &
1506  *( 0.5d0*( 1.0d0+eeta_di(ip, 2)*eeta_lx ) )
1507 
1508  end do
1509 
1510  xxi_lx = xi_lx /dsqrt( 1.0d0/3.0d0 )
1511  eeta_lx = eta_lx/dsqrt( 1.0d0/3.0d0 )
1512 
1513  do ip = 1, npoints_tying(3)
1514 
1515  h(ip, 3) &
1516  = ( 0.5d0*( 1.0d0+xxi_di(ip, 1)*xxi_lx ) ) &
1517  *( 0.5d0*( 1.0d0+eeta_di(ip, 1)*eeta_lx ) )
1518 
1519  end do
1520 
1521  do jsize = 1, ndof*nn
1522 
1523  b(1, jsize) = 0.0d0
1524  b(2, jsize) = 0.0d0
1525  b(3, jsize) = 0.0d0
1526  b(4, jsize) = 0.0d0
1527  b(5, jsize) = 0.0d0
1528 
1529  do ip = 1, npoints_tying(1)
1530 
1531  ! B_as(1, jsize)
1532  b(1, jsize) &
1533  = b(1, jsize)+h(ip, 1)*b_di(1, jsize, ip, 1, ly)
1534  ! B_as(5, jsize)
1535  b(5, jsize) &
1536  = b(5, jsize)+h(ip, 1)*b_di(5, jsize, ip, 1, ly)
1537 
1538  end do
1539 
1540  do ip = 1, npoints_tying(2)
1541 
1542  ! B_as(2, jsize)
1543  b(2, jsize) &
1544  = b(2, jsize)+h(ip, 2)*b_di(2, jsize, ip, 2, ly)
1545  ! B_as(4, jsize)
1546  b(4, jsize) &
1547  = b(4, jsize)+h(ip, 2)*b_di(4, jsize, ip, 2, ly)
1548 
1549  end do
1550 
1551  do ip = 1, npoints_tying(3)
1552 
1553  ! B_as(3, jsize)
1554  b(3, jsize) &
1555  = b(3, jsize)+h(ip, 3)*b_di(3, jsize, ip, 3, ly)
1556 
1557  end do
1558 
1559  end do
1560 
1561  ! MITC3
1562  else if( etype .EQ. fe_mitc3_shell ) then
1563 
1564  do jsize = 1, ndof*nn
1565 
1566  b(4, jsize) = 0.0d0
1567  b(5, jsize) = 0.0d0
1568 
1569  ! B_as(4, jsize)
1570  b(4, jsize) &
1571  = ( 1.0d0-xi_lx )*b_di(4, jsize, 2, 1, ly) &
1572  +xi_lx *b_di(5, jsize, 1, 1, ly) &
1573  +xi_lx *( b_di(4, jsize, 3, 1, ly) &
1574  -b_di(5, jsize, 3, 1, ly) )
1575 
1576  ! B_as(5, jsize)
1577  b(5, jsize) &
1578  = eta_lx*b_di(4, jsize, 2, 1, ly) &
1579  +( 1.0d0-eta_lx )*b_di(5, jsize, 1, 1, ly) &
1580  -eta_lx*( b_di(4, jsize, 3, 1, ly) &
1581  -b_di(5, jsize, 3, 1, ly) )
1582 
1583  end do
1584 
1585  end if
1586 
1587  !--------------------------------------------------
1588 
1589  w_w_w_det = w_w_lx*w_ly*det_weight
1590 
1591  !--------------------------------------------------
1592 
1593  if( present( qf_stress ) .or. add_geo_stiff ) then
1594  if( ishell > 0 ) then
1595  sv_force(1) = element%shell_layer_gausses(ishell)%stress(1)
1596  sv_force(2) = element%shell_layer_gausses(ishell)%stress(2)
1597  sv_force(3) = element%shell_layer_gausses(ishell)%stress(4)
1598  sv_force(4) = element%shell_layer_gausses(ishell)%stress(5)
1599  sv_force(5) = element%shell_layer_gausses(ishell)%stress(6)
1600  else
1601  sv_force(1) = gausses(lx)%stress(1)
1602  sv_force(2) = gausses(lx)%stress(2)
1603  sv_force(3) = gausses(lx)%stress(4)
1604  sv_force(4) = gausses(lx)%stress(5)
1605  sv_force(5) = gausses(lx)%stress(6)
1606  endif
1607  endif
1608 
1609  if( present( qf_stress ) ) then
1610  qf_tmp(1:ndof*nn) = qf_tmp(1:ndof*nn) &
1611  +w_w_w_det*gausses(1)%pMaterial%shell_var(n_layer)%weight &
1612  *matmul( sv_force(1:5), b(1:5, 1:ndof*nn) )
1613  endif
1614 
1615  !--------------------------------------------------
1616 
1617  db(1:5, 1:ndof*nn) = matmul( d, b(1:5, 1:ndof*nn ) )
1618 
1619  if( flag == updatelag .and. etype == fe_mitc4_shell .and. nn == 4 ) then
1620  if( ishell > 0 ) then
1621  stress_old_vec(1:6) = element%shell_layer_gausses(ishell)%stress_bak(1:6)
1622  trace_coeff = shellplanestresstracecoeff(element%shell_layer_gausses(ishell), n_layer)
1623  else
1624  stress_old_vec(1:6) = gausses(lx)%stress_bak(1:6)
1625  trace_coeff = shellplanestresstracecoeff(gausses(lx), n_layer)
1626  endif
1627  call shelladdulobjectivetangent( stress_old_vec, ndof*nn, b, db, trace_coeff=trace_coeff )
1628  endif
1629 
1630  !--------------------------------------------------
1631 
1632  do jsize=1,ndof*nn
1633  do isize=1,ndof*nn
1634  tmpstiff(isize, jsize) &
1635  = tmpstiff(isize, jsize) &
1636  +w_w_w_det*gausses(1)%pMaterial%shell_var(n_layer)%weight*dot_product( b(:, isize), db(:, jsize) )
1637  end do
1638  end do
1639 
1640  !--------------------------------------------------
1641 
1642  ! [ B_{i} ] matrix
1643  b1(:, :) = 0.0d0
1644  b2(:, :) = 0.0d0
1645  b3(:, :) = 0.0d0
1646  cv_w_second(:, :, :) = 0.0d0
1647  do nb = 1, nn
1648 
1649  jsize1 = ndof*(nb-1)+1
1650  jsize2 = ndof*(nb-1)+2
1651  jsize3 = ndof*(nb-1)+3
1652  jsize4 = ndof*(nb-1)+4
1653  jsize5 = ndof*(nb-1)+5
1654  jsize6 = ndof*(nb-1)+6
1655 
1656  b1(1, jsize1) = shapederiv(nb, 1)
1657  b1(2, jsize1) = 0.0d0
1658  b1(3, jsize1) = 0.0d0
1659  b1(1, jsize2) = 0.0d0
1660  b1(2, jsize2) = shapederiv(nb, 1)
1661  b1(3, jsize2) = 0.0d0
1662  b1(1, jsize3) = 0.0d0
1663  b1(2, jsize3) = 0.0d0
1664  b1(3, jsize3) = shapederiv(nb, 1)
1665  b1(1, jsize4) = 0.0d0
1666  b1(2, jsize4) = -dudxi_rot(3, nb)
1667  b1(3, jsize4) = dudxi_rot(2, nb)
1668  b1(1, jsize5) = dudxi_rot(3, nb)
1669  b1(2, jsize5) = 0.0d0
1670  b1(3, jsize5) = -dudxi_rot(1, nb)
1671  b1(1, jsize6) = -dudxi_rot(2, nb)
1672  b1(2, jsize6) = dudxi_rot(1, nb)
1673  b1(3, jsize6) = 0.0d0
1674 
1675  b2(1, jsize1) = shapederiv(nb, 2)
1676  b2(2, jsize1) = 0.0d0
1677  b2(3, jsize1) = 0.0d0
1678  b2(1, jsize2) = 0.0d0
1679  b2(2, jsize2) = shapederiv(nb, 2)
1680  b2(3, jsize2) = 0.0d0
1681  b2(1, jsize3) = 0.0d0
1682  b2(2, jsize3) = 0.0d0
1683  b2(3, jsize3) = shapederiv(nb, 2)
1684  b2(1, jsize4) = 0.0d0
1685  b2(2, jsize4) = -dudeta_rot(3, nb)
1686  b2(3, jsize4) = dudeta_rot(2, nb)
1687  b2(1, jsize5) = dudeta_rot(3, nb)
1688  b2(2, jsize5) = 0.0d0
1689  b2(3, jsize5) = -dudeta_rot(1, nb)
1690  b2(1, jsize6) = -dudeta_rot(2, nb)
1691  b2(2, jsize6) = dudeta_rot(1, nb)
1692  b2(3, jsize6) = 0.0d0
1693 
1694  b3(1, jsize1) = 0.0d0
1695  b3(2, jsize1) = 0.0d0
1696  b3(3, jsize1) = 0.0d0
1697  b3(1, jsize2) = 0.0d0
1698  b3(2, jsize2) = 0.0d0
1699  b3(3, jsize2) = 0.0d0
1700  b3(1, jsize3) = 0.0d0
1701  b3(2, jsize3) = 0.0d0
1702  b3(3, jsize3) = 0.0d0
1703  b3(1, jsize4) = 0.0d0
1704  b3(2, jsize4) = -dudzeta_rot(3, nb)
1705  b3(3, jsize4) = dudzeta_rot(2, nb)
1706  b3(1, jsize5) = dudzeta_rot(3, nb)
1707  b3(2, jsize5) = 0.0d0
1708  b3(3, jsize5) = -dudzeta_rot(1, nb)
1709  b3(1, jsize6) = -dudzeta_rot(2, nb)
1710  b3(2, jsize6) = dudzeta_rot(1, nb)
1711  b3(3, jsize6) = 0.0d0
1712 
1713  if( use_director_tangent ) then
1714  dudxi_rot_deriv(1:3, 1:3) = shapederiv(nb, 1)*zeta_ly &
1715  *a_over_2_v3_deriv(1:3, 1:3, nb)
1716  dudeta_rot_deriv(1:3, 1:3) = shapederiv(nb, 2)*zeta_ly &
1717  *a_over_2_v3_deriv(1:3, 1:3, nb)
1718  dudzeta_rot_deriv(1:3, 1:3) = shapefunc(nb) &
1719  *a_over_2_v3_deriv(1:3, 1:3, nb)
1720  dudxi_rot_second(1:3, 1:3, 1:3) = shapederiv(nb, 1)*zeta_ly &
1721  *a_over_2_v3_second(1:3, 1:3, 1:3, nb)
1722  dudeta_rot_second(1:3, 1:3, 1:3) = shapederiv(nb, 2)*zeta_ly &
1723  *a_over_2_v3_second(1:3, 1:3, 1:3, nb)
1724  dudzeta_rot_second(1:3, 1:3, 1:3) = shapefunc(nb) &
1725  *a_over_2_v3_second(1:3, 1:3, 1:3, nb)
1726 
1727  do m = 1, 3
1728  jsize = ndof*(nb-1)+3+m
1729  b1(1:3, jsize) = dudxi_rot_deriv(1:3, m)
1730  b2(1:3, jsize) = dudeta_rot_deriv(1:3, m)
1731  b3(1:3, jsize) = dudzeta_rot_deriv(1:3, m)
1732  end do
1733  do n = 1, 3
1734  do m = 1, 3
1735  cv12_2 = ( cg1_mat(1)*dudxi_rot_second(2, m, n) &
1736  +cg2_mat(1)*dudeta_rot_second(2, m, n) &
1737  +cg3_mat(1)*dudzeta_rot_second(2, m, n) ) &
1738  -( cg1_mat(2)*dudxi_rot_second(1, m, n) &
1739  +cg2_mat(2)*dudeta_rot_second(1, m, n) &
1740  +cg3_mat(2)*dudzeta_rot_second(1, m, n) )
1741  cv13_2 = ( cg1_mat(1)*dudxi_rot_second(3, m, n) &
1742  +cg2_mat(1)*dudeta_rot_second(3, m, n) &
1743  +cg3_mat(1)*dudzeta_rot_second(3, m, n) ) &
1744  -( cg1_mat(3)*dudxi_rot_second(1, m, n) &
1745  +cg2_mat(3)*dudeta_rot_second(1, m, n) &
1746  +cg3_mat(3)*dudzeta_rot_second(1, m, n) )
1747  cv21_2 = ( cg1_mat(2)*dudxi_rot_second(1, m, n) &
1748  +cg2_mat(2)*dudeta_rot_second(1, m, n) &
1749  +cg3_mat(2)*dudzeta_rot_second(1, m, n) ) &
1750  -( cg1_mat(1)*dudxi_rot_second(2, m, n) &
1751  +cg2_mat(1)*dudeta_rot_second(2, m, n) &
1752  +cg3_mat(1)*dudzeta_rot_second(2, m, n) )
1753  cv23_2 = ( cg1_mat(2)*dudxi_rot_second(3, m, n) &
1754  +cg2_mat(2)*dudeta_rot_second(3, m, n) &
1755  +cg3_mat(2)*dudzeta_rot_second(3, m, n) ) &
1756  -( cg1_mat(3)*dudxi_rot_second(2, m, n) &
1757  +cg2_mat(3)*dudeta_rot_second(2, m, n) &
1758  +cg3_mat(3)*dudzeta_rot_second(2, m, n) )
1759  cv31_2 = ( cg1_mat(3)*dudxi_rot_second(1, m, n) &
1760  +cg2_mat(3)*dudeta_rot_second(1, m, n) &
1761  +cg3_mat(3)*dudzeta_rot_second(1, m, n) ) &
1762  -( cg1_mat(1)*dudxi_rot_second(3, m, n) &
1763  +cg2_mat(1)*dudeta_rot_second(3, m, n) &
1764  +cg3_mat(1)*dudzeta_rot_second(3, m, n) )
1765  cv32_2 = ( cg1_mat(3)*dudxi_rot_second(2, m, n) &
1766  +cg2_mat(3)*dudeta_rot_second(2, m, n) &
1767  +cg3_mat(3)*dudzeta_rot_second(2, m, n) ) &
1768  -( cg1_mat(2)*dudxi_rot_second(3, m, n) &
1769  +cg2_mat(2)*dudeta_rot_second(3, m, n) &
1770  +cg3_mat(2)*dudzeta_rot_second(3, m, n) )
1771  cv_w_second(m, n, nb) = v1_i(1)*cv12_2*v2_i(2) &
1772  +v1_i(1)*cv13_2*v2_i(3) &
1773  +v1_i(2)*cv21_2*v2_i(1) &
1774  +v1_i(2)*cv23_2*v2_i(3) &
1775  +v1_i(3)*cv31_2*v2_i(1) &
1776  +v1_i(3)*cv32_2*v2_i(2)
1777  end do
1778  end do
1779  endif
1780 
1781  end do
1782 
1783  !--------------------------------------------------
1784 
1785  if( flag == updatelag ) then
1786  do jsize = 1, ndof*nn
1787  vol_deriv(jsize) = dot_product(cg1(1:3), b1(1:3, jsize)) &
1788  +dot_product(cg2(1:3), b2(1:3, jsize)) &
1789  +dot_product(cg3(1:3), b3(1:3, jsize))
1790  qf_stress_integrand(jsize) = dot_product(sv_force(1:5), b(1:5, jsize))
1791  end do
1792  endif
1793 
1794  !--------------------------------------------------
1795 
1796  if( add_geo_stiff ) then
1797  if( etype == fe_mitc4_shell .and. ( use_tl_green .or. flag == updatelag ) ) then
1798  ! Match the MITC4 residual integration points.
1799  do jsize=1,ndof*nn
1800  do isize=1,ndof*nn
1801  geo_term = sv_force(1)*dot_product(b1(1:3, isize), b1(1:3, jsize)) &
1802  +sv_force(2)*dot_product(b2(1:3, isize), b2(1:3, jsize)) &
1803  +sv_force(3)*(dot_product(b1(1:3, isize), b2(1:3, jsize)) &
1804  +dot_product(b2(1:3, isize), b1(1:3, jsize))) &
1805  +sv_force(4)*(0.5d0*(1.0d0-xi_lx) &
1806  *(dot_product(bg2_di(1:3, isize, 4, 1, ly), bg3_di(1:3, jsize, 4, 1, ly)) &
1807  +dot_product(bg3_di(1:3, isize, 4, 1, ly), bg2_di(1:3, jsize, 4, 1, ly))) &
1808  +0.5d0*(1.0d0+xi_lx) &
1809  *(dot_product(bg2_di(1:3, isize, 2, 1, ly), bg3_di(1:3, jsize, 2, 1, ly)) &
1810  +dot_product(bg3_di(1:3, isize, 2, 1, ly), bg2_di(1:3, jsize, 2, 1, ly)))) &
1811  +sv_force(5)*(0.5d0*(1.0d0-eta_lx) &
1812  *(dot_product(bg3_di(1:3, isize, 1, 1, ly), bg1_di(1:3, jsize, 1, 1, ly)) &
1813  +dot_product(bg1_di(1:3, isize, 1, 1, ly), bg3_di(1:3, jsize, 1, 1, ly))) &
1814  +0.5d0*(1.0d0+eta_lx) &
1815  *(dot_product(bg3_di(1:3, isize, 3, 1, ly), bg1_di(1:3, jsize, 3, 1, ly)) &
1816  +dot_product(bg1_di(1:3, isize, 3, 1, ly), bg3_di(1:3, jsize, 3, 1, ly))))
1817  tmpstiff(isize, jsize) = tmpstiff(isize, jsize) &
1818  +w_w_w_det*gausses(1)%pMaterial%shell_var(n_layer)%weight &
1819  *geo_term
1820  end do
1821  end do
1822  else
1823  s_global(:, :) = 0.0d0
1824  s_global(:, :) = s_global(:, :) &
1825  +sv_force(1)*outer_product3(e1_hat_mat, e1_hat_mat) &
1826  +sv_force(2)*outer_product3(e2_hat_mat, e2_hat_mat) &
1827  +sv_force(3)*(outer_product3(e1_hat_mat, e2_hat_mat) &
1828  +outer_product3(e2_hat_mat, e1_hat_mat)) &
1829  +sv_force(4)*(outer_product3(e2_hat_mat, e3_hat_mat) &
1830  +outer_product3(e3_hat_mat, e2_hat_mat)) &
1831  +sv_force(5)*(outer_product3(e1_hat_mat, e3_hat_mat) &
1832  +outer_product3(e3_hat_mat, e1_hat_mat))
1833 
1834  bn(1:9, 1:ndof*nn) = 0.0d0
1835  do jsize = 1, ndof*nn
1836  bn(1, jsize) = b1(1, jsize)*cg1_mat(1) + b2(1, jsize)*cg2_mat(1) &
1837  +b3(1, jsize)*cg3_mat(1)
1838  bn(2, jsize) = b1(2, jsize)*cg1_mat(1) + b2(2, jsize)*cg2_mat(1) &
1839  +b3(2, jsize)*cg3_mat(1)
1840  bn(3, jsize) = b1(3, jsize)*cg1_mat(1) + b2(3, jsize)*cg2_mat(1) &
1841  +b3(3, jsize)*cg3_mat(1)
1842  bn(4, jsize) = b1(1, jsize)*cg1_mat(2) + b2(1, jsize)*cg2_mat(2) &
1843  +b3(1, jsize)*cg3_mat(2)
1844  bn(5, jsize) = b1(2, jsize)*cg1_mat(2) + b2(2, jsize)*cg2_mat(2) &
1845  +b3(2, jsize)*cg3_mat(2)
1846  bn(6, jsize) = b1(3, jsize)*cg1_mat(2) + b2(3, jsize)*cg2_mat(2) &
1847  +b3(3, jsize)*cg3_mat(2)
1848  bn(7, jsize) = b1(1, jsize)*cg1_mat(3) + b2(1, jsize)*cg2_mat(3) &
1849  +b3(1, jsize)*cg3_mat(3)
1850  bn(8, jsize) = b1(2, jsize)*cg1_mat(3) + b2(2, jsize)*cg2_mat(3) &
1851  +b3(2, jsize)*cg3_mat(3)
1852  bn(9, jsize) = b1(3, jsize)*cg1_mat(3) + b2(3, jsize)*cg2_mat(3) &
1853  +b3(3, jsize)*cg3_mat(3)
1854  end do
1855 
1856  smat(:, :) = 0.0d0
1857  do j = 1, 3
1858  smat(j , j ) = s_global(1, 1)
1859  smat(j , j+3) = s_global(1, 2)
1860  smat(j , j+6) = s_global(1, 3)
1861  smat(j+3, j ) = s_global(2, 1)
1862  smat(j+3, j+3) = s_global(2, 2)
1863  smat(j+3, j+6) = s_global(2, 3)
1864  smat(j+6, j ) = s_global(3, 1)
1865  smat(j+6, j+3) = s_global(3, 2)
1866  smat(j+6, j+6) = s_global(3, 3)
1867  end do
1868 
1869  sbn(1:9, 1:ndof*nn) = matmul(smat(1:9, 1:9), bn(1:9, 1:ndof*nn))
1870  do jsize=1,ndof*nn
1871  do isize=1,ndof*nn
1872  tmpstiff(isize, jsize) = tmpstiff(isize, jsize) &
1873  +w_w_w_det*gausses(1)%pMaterial%shell_var(n_layer)%weight &
1874  *dot_product(bn(:, isize), sbn(:, jsize))
1875  end do
1876  end do
1877  endif
1878  if( flag == updatelag ) then
1879  do jsize=1,ndof*nn
1880  do isize=1,ndof*nn
1881  tmpstiff(isize, jsize) = tmpstiff(isize, jsize) &
1882  +w_w_w_det*gausses(1)%pMaterial%shell_var(n_layer)%weight &
1883  *qf_stress_integrand(isize)*vol_deriv(jsize)
1884  end do
1885  end do
1886  endif
1887  if( use_director_tangent .and. use_tl_green ) then
1888  ! Director stress stiffness in the director/drilling split.
1889  do nb = 1, nn
1890  drill_axis(1:3) = a_over_2_v3(1:3, nb)
1891  axis_norm = sqrt(sum(drill_axis(1:3)*drill_axis(1:3)))
1892  if( axis_norm > 0.0d0 ) then
1893  drill_axis(1:3) = drill_axis(1:3)/axis_norm
1894  else
1895  drill_axis(1:3) = (/ 0.0d0, 0.0d0, 1.0d0 /)
1896  endif
1897  rot_projector(:, :) = 0.0d0
1898  do i = 1, 3
1899  rot_projector(i, i) = 1.0d0
1900  end do
1901  do i = 1, 3
1902  do j = 1, 3
1903  rot_projector(i, j) = rot_projector(i, j)-drill_axis(i)*drill_axis(j)
1904  end do
1905  end do
1906  do n = 1, 3
1907  do m = 1, 3
1908  hrot(m, n) = dot_product(sv_force(1:5), b2rot(1:5, m, n, nb))
1909  end do
1910  end do
1911  do n = 1, 3
1912  jsize = ndof*(nb-1)+3+n
1913  do m = 1, 3
1914  isize = ndof*(nb-1)+3+m
1915  hess_coeff = 0.0d0
1916  do j = 1, 3
1917  do i = 1, 3
1918  hess_coeff = hess_coeff &
1919  +rot_projector(i, m)*hrot(i, j)*rot_projector(j, n)
1920  end do
1921  end do
1922  drill_coeff = 0.0d0
1923  do i = 1, 3
1924  drill_coeff = drill_coeff+drill_axis(i)*hrot(i, n)
1925  end do
1926  do j = 1, 3
1927  do i = 1, 3
1928  drill_coeff = drill_coeff &
1929  +rot_projector(i, n)*hrot(i, j)*drill_axis(j)
1930  end do
1931  end do
1932  tmpstiff(isize, jsize) = tmpstiff(isize, jsize) &
1933  +w_w_w_det*gausses(1)%pMaterial%shell_var(n_layer)%weight &
1934  *(hess_coeff+drill_axis(m)*drill_coeff)
1935  end do
1936  end do
1937  end do
1938  endif
1939  if( use_director_tangent .and. flag == updatelag ) then
1940  ! Current-frame counterpart of the director/drilling split.
1941  do nb = 1, nn
1942  drill_axis(1:3) = a_over_2_v3(1:3, nb)
1943  axis_norm = sqrt(sum(drill_axis(1:3)*drill_axis(1:3)))
1944  if( axis_norm > 0.0d0 ) then
1945  drill_axis(1:3) = drill_axis(1:3)/axis_norm
1946  else
1947  drill_axis(1:3) = (/ 0.0d0, 0.0d0, 1.0d0 /)
1948  endif
1949  rot_projector(:, :) = 0.0d0
1950  do i = 1, 3
1951  rot_projector(i, i) = 1.0d0
1952  end do
1953  do i = 1, 3
1954  do j = 1, 3
1955  rot_projector(i, j) = rot_projector(i, j)-drill_axis(i)*drill_axis(j)
1956  end do
1957  end do
1958  do n = 1, 3
1959  do m = 1, 3
1960  hrot(m, n) = dot_product(sv_force(1:5), b2rot(1:5, m, n, nb))
1961  end do
1962  end do
1963  do n = 1, 3
1964  jsize = ndof*(nb-1)+3+n
1965  do m = 1, 3
1966  isize = ndof*(nb-1)+3+m
1967  hess_coeff = 0.0d0
1968  do j = 1, 3
1969  do i = 1, 3
1970  hess_coeff = hess_coeff &
1971  +rot_projector(i, m)*hrot(i, j)*rot_projector(j, n)
1972  end do
1973  end do
1974  drill_coeff = 0.0d0
1975  do i = 1, 3
1976  drill_coeff = drill_coeff+drill_axis(i)*hrot(i, n)
1977  end do
1978  do j = 1, 3
1979  do i = 1, 3
1980  drill_coeff = drill_coeff &
1981  +rot_projector(i, n)*hrot(i, j)*drill_axis(j)
1982  end do
1983  end do
1984  tmpstiff(isize, jsize) = tmpstiff(isize, jsize) &
1985  +w_w_w_det*gausses(1)%pMaterial%shell_var(n_layer)%weight &
1986  *(hess_coeff+drill_axis(m)*drill_coeff)
1987  end do
1988  end do
1989  end do
1990  endif
1991  endif
1992 
1993  !--------------------------------------------------
1994 
1995  ! { C_{ij} } vector
1996  do jsize = 1, ndof*nn
1997 
1998  cv12(jsize) = ( cg1_mat(1)*b1(2, jsize) &
1999  +cg2_mat(1)*b2(2, jsize) &
2000  +cg3_mat(1)*b3(2, jsize) ) &
2001  -( cg1_mat(2)*b1(1, jsize) &
2002  +cg2_mat(2)*b2(1, jsize) &
2003  +cg3_mat(2)*b3(1, jsize) )
2004  cv13(jsize) = ( cg1_mat(1)*b1(3, jsize) &
2005  +cg2_mat(1)*b2(3, jsize) &
2006  +cg3_mat(1)*b3(3, jsize) ) &
2007  -( cg1_mat(3)*b1(1, jsize) &
2008  +cg2_mat(3)*b2(1, jsize) &
2009  +cg3_mat(3)*b3(1, jsize) )
2010  cv21(jsize) = ( cg1_mat(2)*b1(1, jsize) &
2011  +cg2_mat(2)*b2(1, jsize) &
2012  +cg3_mat(2)*b3(1, jsize) ) &
2013  -( cg1_mat(1)*b1(2, jsize) &
2014  +cg2_mat(1)*b2(2, jsize) &
2015  +cg3_mat(1)*b3(2, jsize) )
2016  cv23(jsize) = ( cg1_mat(2)*b1(3, jsize) &
2017  +cg2_mat(2)*b2(3, jsize) &
2018  +cg3_mat(2)*b3(3, jsize) ) &
2019  -( cg1_mat(3)*b1(2, jsize) &
2020  +cg2_mat(3)*b2(2, jsize) &
2021  +cg3_mat(3)*b3(2, jsize) )
2022  cv31(jsize) = ( cg1_mat(3)*b1(1, jsize) &
2023  +cg2_mat(3)*b2(1, jsize) &
2024  +cg3_mat(3)*b3(1, jsize) ) &
2025  -( cg1_mat(1)*b1(3, jsize) &
2026  +cg2_mat(1)*b2(3, jsize) &
2027  +cg3_mat(1)*b3(3, jsize) )
2028  cv32(jsize) = ( cg1_mat(3)*b1(2, jsize) &
2029  +cg2_mat(3)*b2(2, jsize) &
2030  +cg3_mat(3)*b3(2, jsize) ) &
2031  -( cg1_mat(2)*b1(3, jsize) &
2032  +cg2_mat(2)*b2(3, jsize) &
2033  +cg3_mat(2)*b3(3, jsize) )
2034 
2035  end do
2036 
2037  !--------------------------------------------------
2038 
2039  ! { Cw } vector
2040  do nb = 1, nn
2041 
2042  do j = 1, ndof
2043 
2044  jsize = ndof*(nb-1)+j
2045 
2046  cv_w(jsize) &
2047  = v1_i(1)*cv12(jsize)*v2_i(2) &
2048  +v1_i(1)*cv13(jsize)*v2_i(3) &
2049  +v1_i(2)*cv21(jsize)*v2_i(1) &
2050  +v1_i(2)*cv23(jsize)*v2_i(3) &
2051  +v1_i(3)*cv31(jsize)*v2_i(1) &
2052  +v1_i(3)*cv32(jsize)*v2_i(2)
2053 
2054  end do
2055 
2056  end do
2057 
2058  ! { Ctheta } vector
2059  do nb = 1, nn
2060 
2061  jsize1 = ndof*(nb-1)+1
2062  jsize2 = ndof*(nb-1)+2
2063  jsize3 = ndof*(nb-1)+3
2064  jsize4 = ndof*(nb-1)+4
2065  jsize5 = ndof*(nb-1)+5
2066  jsize6 = ndof*(nb-1)+6
2067 
2068  cv_theta(jsize1) = 0.0d0
2069  cv_theta(jsize2) = 0.0d0
2070  cv_theta(jsize3) = 0.0d0
2071  if( finite_rotation_director .and. present( nddrill ) ) then
2072  drill_axis(1:3) = a_over_2_v3(1:3, nb)
2073  axis_norm = sqrt(sum(drill_axis(1:3)*drill_axis(1:3)))
2074  if( axis_norm > 0.0d0 ) then
2075  drill_axis(1:3) = drill_axis(1:3)/axis_norm
2076  else
2077  drill_axis(1:3) = v3_i(1:3)
2078  endif
2079  cv_theta(jsize4) = drill_axis(1)*shapefunc(nb)
2080  cv_theta(jsize5) = drill_axis(2)*shapefunc(nb)
2081  cv_theta(jsize6) = drill_axis(3)*shapefunc(nb)
2082  else
2083  cv_theta(jsize4) = v3_i(1)*shapefunc(nb)
2084  cv_theta(jsize5) = v3_i(2)*shapefunc(nb)
2085  cv_theta(jsize6) = v3_i(3)*shapefunc(nb)
2086  endif
2087 
2088  end do
2089 
2090  ! { C } vector
2091  ! drilling stabilization for compatibility rotation
2092  do jsize = 1, ndof*nn
2093 
2094  cv(jsize) = cv_theta(jsize)-0.5d0*cv_w(jsize)
2095 
2096  end do
2097 
2098  !--------------------------------------------------
2099 
2100  cv_disp = dot_product( cv(1:ndof*nn), qf_disp(1:ndof*nn) )
2101  if( finite_rotation_director .and. present( nddrill ) ) then
2102  cv_disp = -0.5d0*dot_product( cv_w(1:ndof*nn), qf_disp(1:ndof*nn) )
2103  do nb = 1, nn
2104  cv_disp = cv_disp + shapefunc(nb)*nddrill(nb)
2105  end do
2106  endif
2107 
2108  if( present( qf_stress ) ) then
2109  qf_tmp(1:ndof*nn) = qf_tmp(1:ndof*nn) &
2110  +w_w_w_det*gausses(1)%pMaterial%shell_var(n_layer)%weight &
2111  *alpha*cv(1:ndof*nn)*cv_disp
2112  endif
2113 
2114  !--------------------------------------------------
2115 
2116  ! [ K L ] matrix
2117  do jsize = 1, ndof*nn
2118  do isize = 1, ndof*nn
2119 
2120  tmpstiff(isize, jsize) &
2121  = tmpstiff(isize, jsize) &
2122  +w_w_w_det*gausses(1)%pMaterial%shell_var(n_layer)%weight*alpha*cv(isize)*cv(jsize)
2123 
2124  end do
2125  end do
2126 
2127  if( use_director_tangent ) then
2128  do nb = 1, nn
2129  do n = 1, 3
2130  jsize = ndof*(nb-1)+3+n
2131  cv_deriv_disp = 0.0d0
2132  do m = 1, 3
2133  isize = ndof*(nb-1)+3+m
2134  cv_deriv = -0.5d0*cv_w_second(m, n, nb)
2135  if( finite_rotation_director .and. present( nddrill ) ) then
2136  drill_axis(1:3) = a_over_2_v3(1:3, nb)
2137  axis_norm = sqrt(sum(drill_axis(1:3)*drill_axis(1:3)))
2138  if( axis_norm > 0.0d0 ) then
2139  drill_axis(1:3) = drill_axis(1:3)/axis_norm
2140  drill_coeff = dot_product(drill_axis(1:3), a_over_2_v3_deriv(1:3, n, nb))
2141  cv_deriv = cv_deriv + shapefunc(nb) &
2142  *(a_over_2_v3_deriv(m, n, nb) - drill_axis(m)*drill_coeff)/axis_norm
2143  endif
2144  endif
2145  cv_deriv_disp = cv_deriv_disp + cv_deriv*qf_disp(isize)
2146  tmpstiff(isize, jsize) = tmpstiff(isize, jsize) &
2147  +w_w_w_det*gausses(1)%pMaterial%shell_var(n_layer)%weight &
2148  *alpha*cv_deriv*cv_disp
2149  end do
2150  do isize = 1, ndof*nn
2151  tmpstiff(isize, jsize) = tmpstiff(isize, jsize) &
2152  +w_w_w_det*gausses(1)%pMaterial%shell_var(n_layer)%weight &
2153  *alpha*cv(isize)*cv_deriv_disp
2154  end do
2155  end do
2156  end do
2157  endif
2158 
2159  !--------------------------------------------------
2160 
2161  end do
2162 
2163  !--------------------------------------------------------
2164 
2165  end do
2166 
2167  !--------------------------------------------------------------
2168 
2169  stiff(1:nn*ndof, 1:nn*ndof) = tmpstiff(1:nn*ndof, 1:nn*ndof)
2170 
2171  !--------------------------------------------------------------------
2172  end do
2173 
2174  ! write(*,"(24E25.16)") stiff
2175 
2176  !******************** Shell-Solid mixed analysis ********************
2177  ! mixglaf = 0 < natural shell (6 dof)
2178  ! mixglaf = 1 < mixed 361 (3*2 dof)*8 nod
2179  ! mixglaf = 2 < mixed 351 (3*2 dof)*6 nod
2180  if( mixflag == 1 )then
2181 
2182  ! write(*,*) 'convert for shell-solid mixed analysis'
2183  sstable(1) = 1
2184  sstable(2) = 2
2185  sstable(3) = 3
2186  sstable(4) = 7
2187  sstable(5) = 8
2188  sstable(6) = 9
2189  sstable(7) = 13
2190  sstable(8) = 14
2191  sstable(9) = 15
2192  sstable(10)= 19
2193  sstable(11)= 20
2194  sstable(12)= 21
2195  sstable(13)= 4
2196  sstable(14)= 5
2197  sstable(15)= 6
2198  sstable(16)= 10
2199  sstable(17)= 11
2200  sstable(18)= 12
2201  sstable(19)= 16
2202  sstable(20)= 17
2203  sstable(21)= 18
2204  sstable(22)= 22
2205  sstable(23)= 23
2206  sstable(24)= 24
2207 
2208  tmpstiff(1:nn*ndof, 1:nn*ndof) = stiff(1:nn*ndof, 1:nn*ndof)
2209 
2210  do i = 1, nn*ndof
2211  do j = 1, nn*ndof
2212  stiff(i,j) = tmpstiff(sstable(i),sstable(j))
2213  enddo
2214  enddo
2215 
2216  if( present( qf_stress ) ) then
2217  qf_mix(1:nn*ndof) = qf_tmp(1:nn*ndof)
2218  do i = 1, nn*ndof
2219  qf_stress(i) = qf_mix(sstable(i))
2220  enddo
2221  endif
2222 
2223  elseif( mixflag == 2 )then
2224  ! write(*,*) 'convert for shell-solid mixed analysis 351'
2225  sstable(1) = 1
2226  sstable(2) = 2
2227  sstable(3) = 3
2228  sstable(4) = 7
2229  sstable(5) = 8
2230  sstable(6) = 9
2231  sstable(7) = 13
2232  sstable(8) = 14
2233  sstable(9) = 15
2234  sstable(10)= 4
2235  sstable(11)= 5
2236  sstable(12)= 6
2237  sstable(13)= 10
2238  sstable(14)= 11
2239  sstable(15)= 12
2240  sstable(16)= 16
2241  sstable(17)= 17
2242  sstable(18)= 18
2243 
2244  tmpstiff(1:nn*ndof, 1:nn*ndof) = stiff(1:nn*ndof, 1:nn*ndof)
2245 
2246  do i = 1, nn*ndof
2247  do j = 1, nn*ndof
2248  stiff(i,j) = tmpstiff(sstable(i),sstable(j))
2249  enddo
2250  enddo
2251 
2252  if( present( qf_stress ) ) then
2253  qf_mix(1:nn*ndof) = qf_tmp(1:nn*ndof)
2254  do i = 1, nn*ndof
2255  qf_stress(i) = qf_mix(sstable(i))
2256  enddo
2257  endif
2258 
2259  else
2260  if( present( qf_stress ) ) qf_stress(1:nn*ndof) = qf_tmp(1:nn*ndof)
2261 
2262  endif
2263 
2264  return
2265 
2266  !####################################################################
2267  end subroutine stf_shell_mitc
2268  !####################################################################
2269 
2270 
2271  !####################################################################
2272  subroutine shelldirectorincrement( theta, director_ref, director_inc )
2273  !####################################################################
2274 
2275  real(kind = kreal), intent(in) :: theta(3)
2276  real(kind = kreal), intent(in) :: director_ref(3)
2277  real(kind = kreal), intent(out) :: director_inc(3)
2278 
2279  real(kind = kreal) :: theta_norm, theta_norm2
2280  real(kind = kreal) :: cross1(3), cross2(3)
2281  real(kind = kreal) :: sin_over_theta, one_minus_cos_over_theta2
2282 
2283  cross1(1) = theta(2)*director_ref(3) - theta(3)*director_ref(2)
2284  cross1(2) = theta(3)*director_ref(1) - theta(1)*director_ref(3)
2285  cross1(3) = theta(1)*director_ref(2) - theta(2)*director_ref(1)
2286 
2287  cross2(1) = theta(2)*cross1(3) - theta(3)*cross1(2)
2288  cross2(2) = theta(3)*cross1(1) - theta(1)*cross1(3)
2289  cross2(3) = theta(1)*cross1(2) - theta(2)*cross1(1)
2290 
2291  theta_norm2 = dot_product( theta(1:3), theta(1:3) )
2292  theta_norm = dsqrt( theta_norm2 )
2293 
2294  if( theta_norm < 1.0d-12 ) then
2295  director_inc(1:3) = cross1(1:3) + 0.5d0*cross2(1:3)
2296  else
2297  sin_over_theta = dsin( theta_norm )/theta_norm
2298  one_minus_cos_over_theta2 = ( 1.0d0-dcos( theta_norm ) )/theta_norm2
2299  director_inc(1:3) = sin_over_theta*cross1(1:3) &
2300  +one_minus_cos_over_theta2*cross2(1:3)
2301  endif
2302 
2303  return
2304 
2305  !####################################################################
2306  end subroutine shelldirectorincrement
2307  !####################################################################
2308 
2309  !####################################################################
2310  pure subroutine shellrotationvectortomatrix( theta, rotmat )
2311  !####################################################################
2312 
2313  real(kind = kreal), intent(in) :: theta(3)
2314  real(kind = kreal), intent(out) :: rotmat(3, 3)
2315 
2316  real(kind = kreal) :: theta_norm, theta_norm2
2317  real(kind = kreal) :: sin_over_theta, one_minus_cos_over_theta2
2318  real(kind = kreal) :: skew(3, 3), skew2(3, 3)
2319  integer :: i, j
2320 
2321  theta_norm2 = dot_product( theta(1:3), theta(1:3) )
2322  theta_norm = dsqrt( theta_norm2 )
2323 
2324  skew(:, :) = 0.0d0
2325  skew(1, 2) = -theta(3)
2326  skew(1, 3) = theta(2)
2327  skew(2, 1) = theta(3)
2328  skew(2, 3) = -theta(1)
2329  skew(3, 1) = -theta(2)
2330  skew(3, 2) = theta(1)
2331 
2332  skew2 = matmul( skew, skew )
2333 
2334  if( theta_norm < 1.0d-12 ) then
2335  sin_over_theta = 1.0d0 - theta_norm2/6.0d0 + theta_norm2*theta_norm2/120.0d0
2336  one_minus_cos_over_theta2 = 0.5d0 - theta_norm2/24.0d0 + theta_norm2*theta_norm2/720.0d0
2337  else
2338  sin_over_theta = dsin( theta_norm )/theta_norm
2339  one_minus_cos_over_theta2 = ( 1.0d0-dcos( theta_norm ) )/theta_norm2
2340  endif
2341 
2342  rotmat(:, :) = sin_over_theta*skew(:, :) + one_minus_cos_over_theta2*skew2(:, :)
2343  do i = 1, 3
2344  rotmat(i, i) = rotmat(i, i) + 1.0d0
2345  end do
2346 
2347  return
2348 
2349  !####################################################################
2350  end subroutine shellrotationvectortomatrix
2351  !####################################################################
2352 
2353 
2354  !####################################################################
2355  pure subroutine shellrotationmatrixtovector( rotmat, theta )
2356  !####################################################################
2357 
2358  real(kind = kreal), intent(in) :: rotmat(3, 3)
2359  real(kind = kreal), intent(out) :: theta(3)
2360 
2361  real(kind = kreal) :: pi, trace_r, cos_angle, angle, sin_angle
2362  real(kind = kreal) :: axis(3), axis_abs
2363 
2364  pi = 4.0d0*datan( 1.0d0 )
2365  trace_r = rotmat(1, 1) + rotmat(2, 2) + rotmat(3, 3)
2366  cos_angle = 0.5d0*( trace_r - 1.0d0 )
2367  cos_angle = max( -1.0d0, min( 1.0d0, cos_angle ) )
2368  angle = dacos( cos_angle )
2369 
2370  theta(1) = rotmat(3, 2) - rotmat(2, 3)
2371  theta(2) = rotmat(1, 3) - rotmat(3, 1)
2372  theta(3) = rotmat(2, 1) - rotmat(1, 2)
2373 
2374  if( angle < 1.0d-12 ) then
2375  theta(1:3) = 0.5d0*theta(1:3)
2376  else if( pi-angle < 1.0d-8 ) then
2377  axis(1) = dsqrt( max( 0.0d0, 0.5d0*( rotmat(1, 1)+1.0d0 ) ) )
2378  axis(2) = dsqrt( max( 0.0d0, 0.5d0*( rotmat(2, 2)+1.0d0 ) ) )
2379  axis(3) = dsqrt( max( 0.0d0, 0.5d0*( rotmat(3, 3)+1.0d0 ) ) )
2380 
2381  if( rotmat(2, 1)+rotmat(1, 2) < 0.0d0 ) axis(2) = -axis(2)
2382  if( rotmat(3, 1)+rotmat(1, 3) < 0.0d0 ) axis(3) = -axis(3)
2383  axis_abs = dsqrt( dot_product( axis(1:3), axis(1:3) ) )
2384  if( axis_abs > 1.0d-12 ) then
2385  theta(1:3) = angle*axis(1:3)/axis_abs
2386  else
2387  theta(1:3) = 0.0d0
2388  endif
2389  else
2390  sin_angle = dsin( angle )
2391  theta(1:3) = angle*theta(1:3)/( 2.0d0*sin_angle )
2392  endif
2393 
2394  return
2395 
2396  !####################################################################
2397  end subroutine shellrotationmatrixtovector
2398  !####################################################################
2399 
2400 
2401  !####################################################################
2402  pure subroutine shellcomposerotationvector( theta_old, theta_inc, theta_new )
2403  !####################################################################
2404 
2405  real(kind = kreal), intent(in) :: theta_old(3)
2406  real(kind = kreal), intent(in) :: theta_inc(3)
2407  real(kind = kreal), intent(out) :: theta_new(3)
2408 
2409  real(kind = kreal) :: rot_old(3, 3), rot_inc(3, 3), rot_new(3, 3)
2410 
2411  call shellrotationvectortomatrix( theta_old, rot_old )
2412  call shellrotationvectortomatrix( theta_inc, rot_inc )
2413  rot_new = matmul( rot_inc, rot_old )
2414  call shellrotationmatrixtovector( rot_new, theta_new )
2415 
2416  return
2417 
2418  !####################################################################
2419  end subroutine shellcomposerotationvector
2420  !####################################################################
2421 
2422 
2423  !####################################################################
2424  pure subroutine shellrelativerotationvector( theta_old, theta_target, theta_inc )
2425  !####################################################################
2426 
2427  real(kind = kreal), intent(in) :: theta_old(3)
2428  real(kind = kreal), intent(in) :: theta_target(3)
2429  real(kind = kreal), intent(out) :: theta_inc(3)
2430 
2431  real(kind = kreal) :: rot_old(3, 3), rot_target(3, 3), rot_inc(3, 3)
2432 
2433  call shellrotationvectortomatrix( theta_old, rot_old )
2434  call shellrotationvectortomatrix( theta_target, rot_target )
2435  rot_inc = matmul( rot_target, transpose( rot_old ) )
2436  call shellrotationmatrixtovector( rot_inc, theta_inc )
2437 
2438  return
2439 
2440  !####################################################################
2441  end subroutine shellrelativerotationvector
2442  !####################################################################
2443 
2444 
2445  !####################################################################
2446  pure subroutine shellorthonormalizetriad( triad_in, triad_out )
2447  !####################################################################
2448 
2449  real(kind = kreal), intent(in) :: triad_in(3, 3)
2450  real(kind = kreal), intent(out) :: triad_out(3, 3)
2451 
2452  real(kind = kreal) :: e1(3), e2(3), e3(3), normv
2453 
2454  e1(1:3) = triad_in(1:3, 1)
2455  e3(1:3) = triad_in(1:3, 3)
2456 
2457  normv = dsqrt( dot_product( e3(1:3), e3(1:3) ) )
2458  if( normv < 1.0d-14 ) then
2459  e3(1:3) = (/ 0.0d0, 0.0d0, 1.0d0 /)
2460  else
2461  e3(1:3) = e3(1:3)/normv
2462  endif
2463 
2464  e1(1:3) = e1(1:3) - dot_product( e1(1:3), e3(1:3) )*e3(1:3)
2465  normv = dsqrt( dot_product( e1(1:3), e1(1:3) ) )
2466  if( normv < 1.0d-14 ) then
2467  if( dabs(e3(1)) < 0.9d0 ) then
2468  e1(1:3) = (/ 1.0d0, 0.0d0, 0.0d0 /)
2469  else
2470  e1(1:3) = (/ 0.0d0, 1.0d0, 0.0d0 /)
2471  endif
2472  e1(1:3) = e1(1:3) - dot_product( e1(1:3), e3(1:3) )*e3(1:3)
2473  normv = dsqrt( dot_product( e1(1:3), e1(1:3) ) )
2474  endif
2475  e1(1:3) = e1(1:3)/normv
2476 
2477  e2(1) = e3(2)*e1(3)-e3(3)*e1(2)
2478  e2(2) = e3(3)*e1(1)-e3(1)*e1(3)
2479  e2(3) = e3(1)*e1(2)-e3(2)*e1(1)
2480 
2481  triad_out(1:3, 1) = e1(1:3)
2482  triad_out(1:3, 2) = e2(1:3)
2483  triad_out(1:3, 3) = e3(1:3)
2484 
2485  return
2486 
2487  !####################################################################
2488  end subroutine shellorthonormalizetriad
2489  !####################################################################
2490 
2491 
2492  !####################################################################
2493  pure subroutine shellupdatetriadwithincrement( triad_old, drill_old, theta_inc, triad_new, drill_new )
2494  !####################################################################
2495 
2496  real(kind = kreal), intent(in) :: triad_old(3, 3)
2497  real(kind = kreal), intent(in) :: drill_old
2498  real(kind = kreal), intent(in) :: theta_inc(3)
2499  real(kind = kreal), intent(out) :: triad_new(3, 3)
2500  real(kind = kreal), intent(out) :: drill_new
2501 
2502  real(kind = kreal) :: triad_base(3, 3)
2503  real(kind = kreal) :: director(3), theta_step(3), theta_phys(3)
2504  real(kind = kreal) :: drill_acc, drill_inc, theta_norm
2505  real(kind = kreal) :: rot_inc(3, 3), rot_full(3, 3), triad_trial(3, 3)
2506  integer :: isub, nsub
2507 
2508  call shellorthonormalizetriad( triad_old, triad_base )
2509 
2510  theta_norm = sqrt( sum( theta_inc(1:3)*theta_inc(1:3) ) )
2511  nsub = max( 1, ceiling( theta_norm/5.0d-2 ) )
2512  theta_step(1:3) = theta_inc(1:3)/dble(nsub)
2513  drill_acc = drill_old
2514 
2515  do isub = 1, nsub
2516  director(1:3) = triad_base(1:3, 3)
2517  drill_inc = dot_product( theta_step(1:3), director(1:3) )
2518  theta_phys(1:3) = theta_step(1:3) - drill_inc*director(1:3)
2519 
2520  call shellrotationvectortomatrix( theta_phys, rot_inc )
2521  triad_trial = matmul( rot_inc, triad_base )
2522  call shellorthonormalizetriad( triad_trial, triad_base )
2523  drill_acc = drill_acc + drill_inc
2524  end do
2525 
2526  ! exact director update with separated drilling rotation
2527  call shellrotationvectortomatrix( theta_inc, rot_full )
2528  triad_trial(1:3, 1:3) = triad_base(1:3, 1:3)
2529  triad_trial(1:3, 3) = matmul( rot_full, triad_old(1:3, 3) )
2530  call shellorthonormalizetriad( triad_trial, triad_new )
2531  drill_new = drill_acc
2532 
2533  return
2534 
2535  !####################################################################
2536  end subroutine shellupdatetriadwithincrement
2537  !####################################################################
2538 
2539 
2540  !####################################################################
2541  subroutine shellcomposenodaldisplacement( ndof, nn, disp_old, disp_inc, disp_new )
2542  !####################################################################
2543 
2544  integer(kind = kint), intent(in) :: ndof
2545  integer(kind = kint), intent(in) :: nn
2546  real(kind = kreal), intent(in) :: disp_old(:, :)
2547  real(kind = kreal), intent(in) :: disp_inc(:, :)
2548  real(kind = kreal), intent(out) :: disp_new(6, nn)
2549 
2550  integer :: i, ndof_copy
2551 
2552  disp_new(:, :) = 0.0d0
2553  ndof_copy = min( ndof, 6 )
2554 
2555  do i = 1, nn
2556  disp_new(1:min(3, ndof_copy), i) = disp_old(1:min(3, ndof_copy), i) &
2557  + disp_inc(1:min(3, ndof_copy), i)
2558  if( ndof_copy >= 6 ) then
2559  call shellcomposerotationvector( disp_old(4:6, i), disp_inc(4:6, i), disp_new(4:6, i) )
2560  else if( ndof_copy > 3 ) then
2561  disp_new(4:ndof_copy, i) = disp_old(4:ndof_copy, i) + disp_inc(4:ndof_copy, i)
2562  endif
2563  end do
2564 
2565  return
2566 
2567  !####################################################################
2568  end subroutine shellcomposenodaldisplacement
2569  !####################################################################
2570 
2571 
2572  !####################################################################
2573  pure subroutine shelldirectorincrementalderiv( director_current, director_deriv )
2574  !####################################################################
2575 
2576  real(kind = kreal), intent(in) :: director_current(3)
2577  real(kind = kreal), intent(out) :: director_deriv(3, 3)
2578 
2579  integer :: i
2580  real(kind = kreal) :: basis(3)
2581 
2582  do i = 1, 3
2583  basis(1:3) = 0.0d0
2584  basis(i) = 1.0d0
2585 
2586  director_deriv(1, i) = basis(2)*director_current(3) - basis(3)*director_current(2)
2587  director_deriv(2, i) = basis(3)*director_current(1) - basis(1)*director_current(3)
2588  director_deriv(3, i) = basis(1)*director_current(2) - basis(2)*director_current(1)
2589  end do
2590 
2591  return
2592 
2593  !####################################################################
2594  end subroutine shelldirectorincrementalderiv
2595  !####################################################################
2596 
2597 
2598  !####################################################################
2599  pure subroutine shelldirectorincrementalsecondderiv( director_current, director_second )
2600  !####################################################################
2601 
2602  real(kind = kreal), intent(in) :: director_current(3)
2603  real(kind = kreal), intent(out) :: director_second(3, 3, 3)
2604 
2605  integer :: m, n
2606  real(kind = kreal) :: basis_m(3), basis_n(3)
2607  real(kind = kreal) :: cross_n(3), cross_mn(3), cross_m(3), cross_nm(3)
2608 
2609  do n = 1, 3
2610  basis_n(1:3) = 0.0d0
2611  basis_n(n) = 1.0d0
2612  cross_n(1) = basis_n(2)*director_current(3) - basis_n(3)*director_current(2)
2613  cross_n(2) = basis_n(3)*director_current(1) - basis_n(1)*director_current(3)
2614  cross_n(3) = basis_n(1)*director_current(2) - basis_n(2)*director_current(1)
2615 
2616  do m = 1, 3
2617  basis_m(1:3) = 0.0d0
2618  basis_m(m) = 1.0d0
2619  cross_m(1) = basis_m(2)*director_current(3) - basis_m(3)*director_current(2)
2620  cross_m(2) = basis_m(3)*director_current(1) - basis_m(1)*director_current(3)
2621  cross_m(3) = basis_m(1)*director_current(2) - basis_m(2)*director_current(1)
2622 
2623  cross_mn(1) = basis_m(2)*cross_n(3) - basis_m(3)*cross_n(2)
2624  cross_mn(2) = basis_m(3)*cross_n(1) - basis_m(1)*cross_n(3)
2625  cross_mn(3) = basis_m(1)*cross_n(2) - basis_m(2)*cross_n(1)
2626  cross_nm(1) = basis_n(2)*cross_m(3) - basis_n(3)*cross_m(2)
2627  cross_nm(2) = basis_n(3)*cross_m(1) - basis_n(1)*cross_m(3)
2628  cross_nm(3) = basis_n(1)*cross_m(2) - basis_n(2)*cross_m(1)
2629 
2630  director_second(1:3, m, n) = 0.5d0*( cross_mn(1:3) + cross_nm(1:3) )
2631  end do
2632  end do
2633 
2634  return
2635 
2636  !####################################################################
2637  end subroutine shelldirectorincrementalsecondderiv
2638  !####################################################################
2639 
2640  !####################################################################
2641  subroutine elementstress_shell_mitc &
2642  (etype, nn, ndof, ecoord, gausses, edisp, &
2643  strain, stress, thick, zeta, n_layer, n_totlyr, surface_gauss_points, &
2644  local_strain, local_stress, local_stress_override, nddirector, ndrefdirector, &
2645  ndbase_disp)
2646  !####################################################################
2647 
2648  use m_fstr, only: opsstype, kopss_solution
2649  use m_utilities, only: get_principal
2650  use mmechgauss
2651  use m_matmatrix
2652  use mmaterial, only: infinitesimal, totallag, updatelag, iselastic
2653 
2654  !--------------------------------------------------------------------
2655 
2656  integer(kind = kint), intent(in) :: etype
2657  integer(kind = kint), intent(in) :: nn
2658  integer(kind = kint), intent(in) :: ndof
2659  real(kind = kreal), intent(in) :: ecoord(3, nn)
2660  type(tgaussstatus), intent(in) :: gausses(:)
2661  real(kind = kreal), intent(in) :: edisp(6, nn)
2662  real(kind = kreal), intent(out) :: strain(:,:)
2663  real(kind = kreal), intent(out) :: stress(:,:)
2664  real(kind = kreal), intent(in) :: thick
2665  real(kind = kreal), intent(in) :: zeta
2666  logical, intent(in), optional :: surface_gauss_points
2667  real(kind = kreal), intent(out), optional :: local_strain(:,:)
2668  real(kind = kreal), intent(out), optional :: local_stress(:,:)
2669  real(kind = kreal), intent(in), optional :: local_stress_override(:,:)
2670  real(kind = kreal), intent(in), optional :: nddirector(3, nn)
2671  real(kind = kreal), intent(in), optional :: ndrefdirector(3, nn)
2672  real(kind = kreal), intent(in), optional :: ndbase_disp(6, nn)
2673 
2674  !--------------------------------------------------------------------
2675 
2676  integer :: i
2677  integer :: flag
2678  integer(kind=kint) :: ierr_quad
2679  integer :: lx, npoints
2680  integer :: fetype
2681  integer :: ntying
2682  integer :: npoints_tying(3)
2683  integer :: it, ip
2684  integer :: na, nb
2685  integer :: n_layer, n_totlyr
2686 
2687  real(kind = kreal) :: d(5, 5)
2688  real(kind = kreal) :: elem(3, nn)
2689  real(kind = kreal) :: xi_lx, eta_lx, zeta_ly
2690  real(kind = kreal) :: naturalcoord(2)
2691  real(kind = kreal) :: tpcoord(6, 2, 3)
2692  real(kind = kreal) :: nncoord(nn, 2)
2693  real(kind = kreal) :: shapefunc(nn)
2694  real(kind = kreal) :: shapederiv(nn, 2)
2695  real(kind = kreal) :: alpha
2696  real(kind = kreal) :: xxi_lx, eeta_lx
2697  real(kind = kreal) :: xxi_di(6, 3), eeta_di(6, 3)
2698  real(kind = kreal) :: h(nn, 3)
2699  real(kind = kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
2700  real(kind = kreal) :: v1_abs, v2_abs, v3_abs
2701  real(kind = kreal) :: a_over_2_v3(3, nn)
2702  real(kind = kreal) :: a_over_2_theta_cross_v3(3, nn)
2703  real(kind = kreal) :: u_rot(3, nn)
2704  real(kind = kreal) :: theta(3, nn)
2705  real(kind = kreal) :: dudxi(3), dudeta(3), dudzeta(3)
2706  real(kind = kreal) :: dudxi_trans(3), dudeta_trans(3)
2707  real(kind = kreal) :: g1_cur(3), g2_cur(3), g3_cur(3)
2708  real(kind = kreal) :: g1_ref(3), g2_ref(3), g3_ref(3)
2709  real(kind = kreal) :: cg1_ref(3), cg2_ref(3), cg3_ref(3)
2710  real(kind = kreal) :: det_cur, det_ref, jac
2711  real(kind = kreal) :: dudxi_rot(3, nn), dudeta_rot(3, nn), &
2712  dudzeta_rot(3, nn)
2713  real(kind = kreal) :: g1(3), g2(3), g3(3)
2714  real(kind = kreal) :: g3_abs
2715  real(kind = kreal) :: e_0(3)
2716  real(kind = kreal) :: cg1(3), cg2(3), cg3(3)
2717  real(kind = kreal) :: det
2718  real(kind = kreal) :: det_cg3(3)
2719  real(kind = kreal) :: det_inv
2720  real(kind = kreal) :: det_cg3_abs
2721  real(kind = kreal) :: e1_hat(3), e2_hat(3), e3_hat(3)
2722  real(kind = kreal) :: e1_hat_abs, e2_hat_abs
2723  real(kind = kreal) :: e1_hat_mat(3), e2_hat_mat(3), e3_hat_mat(3)
2724  real(kind = kreal) :: cg1_mat(3), cg2_mat(3), cg3_mat(3)
2725  real(kind = kreal) :: det_mat
2726  real(kind = kreal) :: e11, e22, e12_2, e23_2, e31_2
2727  real(kind = kreal) :: e11_di(6, 3), e22_di(6, 3), &
2728  e12_di_2(6, 3), e23_di_2(6, 3), &
2729  e31_di_2(6, 3)
2730  real(kind = kreal) :: e(3, 3), ev(5)
2731  real(kind = kreal) :: s(3, 3), sv(5)
2732  real(kind = kreal) :: stretch_b(3, 3), tensor(6), eigval(3), princ(3, 3)
2733  real(kind = kreal) :: cauchy(3, 3), logstrain(3, 3), cg_metric(3, 3)
2734  real(kind = kreal) :: eig_norm
2735  logical :: use_surface_gauss
2736  logical :: finite_rotation_director
2737  logical :: use_gl_strain
2738 
2739 
2740  zeta_ly = 0.0d0
2741  use_surface_gauss = .false.
2742  if( present( surface_gauss_points ) ) use_surface_gauss = surface_gauss_points
2743  flag = gausses(1)%pMaterial%nlgeom_flag
2744  finite_rotation_director = ( flag == totallag .or. flag == updatelag ) &
2745  .and. shellsupportsfiniterotationkinematics( etype, nn ) &
2746  .and. iselastic(gausses(1)%pMaterial%mtype)
2747  ! Green-Lagrange strain and spatial output for elastic finite-rotation shell elements.
2748  use_gl_strain = flag == totallag .and. shellsupportsfiniterotationkinematics( etype, nn ) &
2749  .and. iselastic(gausses(1)%pMaterial%mtype)
2750 
2751  !--------------------------------------------------------------------
2752 
2753  ! for lamina stress
2754 
2755  ! MITC4
2756  if( etype .EQ. fe_mitc4_shell ) then
2757 
2758  fetype = fe_mitc4_shell
2759 
2760  ntying = 1
2761  npoints_tying(1)= 4
2762 
2763  ! MITC9
2764  else if( etype .EQ. fe_mitc9_shell ) then
2765 
2766  fetype = fe_mitc9_shell
2767 
2768  ntying = 3
2769  npoints_tying(1)= 6
2770  npoints_tying(2)= 6
2771  npoints_tying(3)= 4
2772 
2773  ! MITC3
2774  else if( etype .EQ. fe_mitc3_shell ) then
2775 
2776  fetype = fe_mitc3_shell
2777 
2778  ntying = 1
2779  npoints_tying(1)= 3
2780 
2781  end if
2782 
2783  !--------------------------------------------------------------------
2784 
2785  elem(:, :) = ecoord(:, :)
2786  if( flag == updatelag .and. present( ndbase_disp ) ) elem(:, :) = elem(:, :) + ndbase_disp(1:3, :)
2787  if( flag == updatelag ) elem(:, :) = elem(:, :) + 0.5d0*edisp(1:3, :)
2788 
2789  !--------------------------------------------------------------------
2790 
2791  do na = 1, nn
2792 
2793  theta(1, na) = edisp(4, na)
2794  theta(2, na) = edisp(5, na)
2795  theta(3, na) = edisp(6, na)
2796 
2797  end do
2798 
2799  !-------------------------------------------------------------------
2800 
2801  ! xi-coordinate at a node in a local element
2802  ! eta-coordinate at a node in a local element
2803  call getnodalnaturalcoord(fetype, nncoord)
2804 
2805  !-------------------------------------------------------------------
2806 
2807  ! MITC4
2808  if( etype .EQ. fe_mitc4_shell ) then
2809 
2810  !--------------------------------------------------------
2811 
2812  ! xi-coordinate at a tying point in a local element
2813  tpcoord(1, 1, 1) = 0.0d0
2814  tpcoord(2, 1, 1) = 1.0d0
2815  tpcoord(3, 1, 1) = 0.0d0
2816  tpcoord(4, 1, 1) = -1.0d0
2817  ! eta-coordinate at a tying point in a local element
2818  tpcoord(1, 2, 1) = -1.0d0
2819  tpcoord(2, 2, 1) = 0.0d0
2820  tpcoord(3, 2, 1) = 1.0d0
2821  tpcoord(4, 2, 1) = 0.0d0
2822 
2823  !--------------------------------------------------------
2824 
2825  ! MITC9
2826  else if( etype .EQ. fe_mitc9_shell ) then
2827 
2828  !--------------------------------------------------------
2829 
2830  ! xi-coordinate at a tying point in a local element
2831  tpcoord(1, 1, 1) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
2832  tpcoord(2, 1, 1) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
2833  tpcoord(3, 1, 1) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
2834  tpcoord(4, 1, 1) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
2835  tpcoord(5, 1, 1) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
2836  tpcoord(6, 1, 1) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
2837  ! eta-coordinate at a tying point in a local element
2838  tpcoord(1, 2, 1) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
2839  tpcoord(2, 2, 1) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
2840  tpcoord(3, 2, 1) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
2841  tpcoord(4, 2, 1) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
2842  tpcoord(5, 2, 1) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
2843  tpcoord(6, 2, 1) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
2844 
2845  ! xi-coordinate at a tying point in a local element
2846  tpcoord(1, 1, 2) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
2847  tpcoord(2, 1, 2) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
2848  tpcoord(3, 1, 2) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
2849  tpcoord(4, 1, 2) = 1.0d0*dsqrt( 3.0d0/5.0d0 )
2850  tpcoord(5, 1, 2) = 0.0d0*dsqrt( 3.0d0/5.0d0 )
2851  tpcoord(6, 1, 2) = -1.0d0*dsqrt( 3.0d0/5.0d0 )
2852  ! eta-coordinate at a tying point in a local element
2853  tpcoord(1, 2, 2) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
2854  tpcoord(2, 2, 2) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
2855  tpcoord(3, 2, 2) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
2856  tpcoord(4, 2, 2) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
2857  tpcoord(5, 2, 2) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
2858  tpcoord(6, 2, 2) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
2859 
2860  ! xi-coordinate at a tying point in a local element
2861  tpcoord(1, 1, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
2862  tpcoord(2, 1, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
2863  tpcoord(3, 1, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
2864  tpcoord(4, 1, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
2865  ! eta-coordinate at a tying point in a local element
2866  tpcoord(1, 2, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
2867  tpcoord(2, 2, 3) = -1.0d0*dsqrt( 1.0d0/3.0d0 )
2868  tpcoord(3, 2, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
2869  tpcoord(4, 2, 3) = 1.0d0*dsqrt( 1.0d0/3.0d0 )
2870 
2871  !--------------------------------------------------------
2872 
2873  ! Xi-coordinate at a tying point in a local element
2874  xxi_di(1, 1) = -1.0d0
2875  xxi_di(2, 1) = 1.0d0
2876  xxi_di(3, 1) = 1.0d0
2877  xxi_di(4, 1) = -1.0d0
2878  xxi_di(5, 1) = 1.0d0
2879  xxi_di(6, 1) = -1.0d0
2880  ! Eta-coordinate at a tying point in a local element
2881  eeta_di(1, 1) = -1.0d0
2882  eeta_di(2, 1) = -1.0d0
2883  eeta_di(3, 1) = 1.0d0
2884  eeta_di(4, 1) = 1.0d0
2885  eeta_di(5, 1) = 0.0d0
2886  eeta_di(6, 1) = 0.0d0
2887 
2888  ! Xi-coordinate at a tying point in a local element
2889  xxi_di(1, 2) = -1.0d0
2890  xxi_di(2, 2) = 0.0d0
2891  xxi_di(3, 2) = 1.0d0
2892  xxi_di(4, 2) = 1.0d0
2893  xxi_di(5, 2) = 0.0d0
2894  xxi_di(6, 2) = -1.0d0
2895  ! Eta-coordinate at a tying point in a local element
2896  eeta_di(1, 2) = -1.0d0
2897  eeta_di(2, 2) = -1.0d0
2898  eeta_di(3, 2) = -1.0d0
2899  eeta_di(4, 2) = 1.0d0
2900  eeta_di(5, 2) = 1.0d0
2901  eeta_di(6, 2) = 1.0d0
2902 
2903  !--------------------------------------------------------
2904 
2905  ! MITC3
2906  else if( etype .EQ. fe_mitc3_shell ) then
2907 
2908  !--------------------------------------------------------
2909 
2910  ! xi-coordinate at a tying point in a local element
2911  tpcoord(1, 1, 1) = 0.5d0
2912  tpcoord(2, 1, 1) = 0.0d0
2913  tpcoord(3, 1, 1) = 0.5d0
2914  ! eta-coordinate at a tying point in a local element
2915  tpcoord(1, 2, 1) = 0.0d0
2916  tpcoord(2, 2, 1) = 0.5d0
2917  tpcoord(3, 2, 1) = 0.5d0
2918 
2919  !--------------------------------------------------------
2920 
2921  end if
2922 
2923  !--------------------------------------------------------------------
2924 
2925  ! xi-coordinate at the center point in a local element
2926  ! eta-coordinate at the center point in a local element
2927  naturalcoord(1) = 0.0d0
2928  naturalcoord(2) = 0.0d0
2929 
2930  call getshapederiv(fetype, naturalcoord, shapederiv)
2931 
2932  !--------------------------------------------------------------
2933 
2934  ! Covariant basis vector
2935  do i = 1, 3
2936 
2937  g1(i) = 0.0d0
2938 
2939  do na = 1, nn
2940 
2941  g1(i) = g1(i)+shapederiv(na, 1) &
2942  *elem(i, na)
2943 
2944  end do
2945 
2946  end do
2947 
2948  e_0(1) = g1(1)
2949  e_0(2) = g1(2)
2950  e_0(3) = g1(3)
2951 
2952  !--------------------------------------------------------------
2953 
2954  do nb = 1, nn
2955 
2956  !--------------------------------------------------------
2957 
2958  naturalcoord(1) = nncoord(nb, 1)
2959  naturalcoord(2) = nncoord(nb, 2)
2960 
2961  call getshapederiv(fetype, naturalcoord, shapederiv)
2962 
2963  !--------------------------------------------------------
2964 
2965  ! Covariant basis vector
2966  do i = 1, 3
2967 
2968  g1(i) = 0.0d0
2969  g2(i) = 0.0d0
2970 
2971  do na = 1, nn
2972 
2973  g1(i) = g1(i)+shapederiv(na, 1) &
2974  *elem(i, na)
2975  g2(i) = g2(i)+shapederiv(na, 2) &
2976  *elem(i, na)
2977 
2978  end do
2979 
2980  end do
2981 
2982  !--------------------------------------------------------
2983 
2984  det_cg3(1) = g1(2)*g2(3)-g1(3)*g2(2)
2985  det_cg3(2) = g1(3)*g2(1)-g1(1)*g2(3)
2986  det_cg3(3) = g1(1)*g2(2)-g1(2)*g2(1)
2987 
2988  det_cg3_abs = dsqrt( det_cg3(1)*det_cg3(1) &
2989  +det_cg3(2)*det_cg3(2) &
2990  +det_cg3(3)*det_cg3(3) )
2991 
2992  v3(1, nb) = det_cg3(1)/det_cg3_abs
2993  v3(2, nb) = det_cg3(2)/det_cg3_abs
2994  v3(3, nb) = det_cg3(3)/det_cg3_abs
2995 
2996  !--------------------------------------------------------
2997 
2998  v2(1, nb) = v3(2, nb)*e_0(3)-v3(3, nb)*e_0(2)
2999  v2(2, nb) = v3(3, nb)*e_0(1)-v3(1, nb)*e_0(3)
3000  v2(3, nb) = v3(1, nb)*e_0(2)-v3(2, nb)*e_0(1)
3001 
3002  v2_abs = dsqrt( v2(1, nb)*v2(1, nb) &
3003  +v2(2, nb)*v2(2, nb) &
3004  +v2(3, nb)*v2(3, nb) )
3005 
3006  if( v2_abs .GT. 1.0d-15 ) then
3007 
3008  v2(1, nb) = v2(1, nb)/v2_abs
3009  v2(2, nb) = v2(2, nb)/v2_abs
3010  v2(3, nb) = v2(3, nb)/v2_abs
3011 
3012  v1(1, nb) = v2(2, nb)*v3(3, nb) &
3013  -v2(3, nb)*v3(2, nb)
3014  v1(2, nb) = v2(3, nb)*v3(1, nb) &
3015  -v2(1, nb)*v3(3, nb)
3016  v1(3, nb) = v2(1, nb)*v3(2, nb) &
3017  -v2(2, nb)*v3(1, nb)
3018 
3019  v1_abs = dsqrt( v1(1, nb)*v1(1, nb) &
3020  +v1(2, nb)*v1(2, nb) &
3021  +v1(3, nb)*v1(3, nb) )
3022 
3023  v1(1, nb) = v1(1, nb)/v1_abs
3024  v1(2, nb) = v1(2, nb)/v1_abs
3025  v1(3, nb) = v1(3, nb)/v1_abs
3026 
3027  else
3028 
3029  v1(1, nb) = 0.0d0
3030  v1(2, nb) = 0.0d0
3031  v1(3, nb) = -1.0d0
3032 
3033  v2(1, nb) = 0.0d0
3034  v2(2, nb) = 1.0d0
3035  v2(3, nb) = 0.0d0
3036 
3037  end if
3038 
3039  !--------------------------------------------------------
3040 
3041  v3(1, nb) = v1(2, nb)*v2(3, nb) &
3042  -v1(3, nb)*v2(2, nb)
3043  v3(2, nb) = v1(3, nb)*v2(1, nb) &
3044  -v1(1, nb)*v2(3, nb)
3045  v3(3, nb) = v1(1, nb)*v2(2, nb) &
3046  -v1(2, nb)*v2(1, nb)
3047 
3048  v3_abs = dsqrt( v3(1, nb)*v3(1, nb) &
3049  +v3(2, nb)*v3(2, nb) &
3050  +v3(3, nb)*v3(3, nb) )
3051 
3052  v3(1, nb) = v3(1, nb)/v3_abs
3053  v3(2, nb) = v3(2, nb)/v3_abs
3054  v3(3, nb) = v3(3, nb)/v3_abs
3055 
3056  !--------------------------------------------------------
3057 
3058  a_over_2_v3(1, nb) = 0.5d0*thick*v3(1, nb)
3059  a_over_2_v3(2, nb) = 0.5d0*thick*v3(2, nb)
3060  a_over_2_v3(3, nb) = 0.5d0*thick*v3(3, nb)
3061  if( finite_rotation_director .and. present( ndrefdirector ) ) then
3062  ! fixed nodal reference director
3063  a_over_2_v3(1:3, nb) = ndrefdirector(1:3, nb)
3064  endif
3065 
3066  !--------------------------------------------------------
3067 
3068  if( finite_rotation_director ) then
3069  if( present( nddirector ) ) then
3070  a_over_2_theta_cross_v3(1:3, nb) = nddirector(1:3, nb) - a_over_2_v3(1:3, nb)
3071  if( flag == updatelag ) then
3072  a_over_2_v3(1:3, nb) = a_over_2_v3(1:3, nb) + 0.5d0*a_over_2_theta_cross_v3(1:3, nb)
3073  else
3074  a_over_2_v3(1:3, nb) = nddirector(1:3, nb)
3075  endif
3076  else
3077  call shelldirectorincrement( theta(1:3, nb), a_over_2_v3(1:3, nb), &
3078  a_over_2_theta_cross_v3(1:3, nb) )
3079  a_over_2_v3(1:3, nb) = a_over_2_v3(1:3, nb) &
3080  + a_over_2_theta_cross_v3(1:3, nb)
3081  endif
3082  else
3083  a_over_2_theta_cross_v3(1, nb) &
3084  = theta(2, nb)*a_over_2_v3(3, nb) &
3085  -theta(3, nb)*a_over_2_v3(2, nb)
3086  a_over_2_theta_cross_v3(2, nb) &
3087  = theta(3, nb)*a_over_2_v3(1, nb) &
3088  -theta(1, nb)*a_over_2_v3(3, nb)
3089  a_over_2_theta_cross_v3(3, nb) &
3090  = theta(1, nb)*a_over_2_v3(2, nb) &
3091  -theta(2, nb)*a_over_2_v3(1, nb)
3092  endif
3093 
3094  !--------------------------------------------------------
3095 
3096  end do
3097 
3098  !--------------------------------------------------------------------
3099  ! Modified stress in laminated shell
3100  !--------------------------------------------------------------------
3101  ! MITC4
3102  if( etype .EQ. fe_mitc4_shell ) then
3103 
3104  zeta_ly = 0.0d0
3105 
3106  ! MITC9
3107  else if( etype .EQ. fe_mitc9_shell ) then
3108 
3109  zeta_ly = zeta
3110 
3111  ! MITC3
3112  else if( etype .EQ. fe_mitc3_shell ) then
3113 
3114  zeta_ly = 0.0d0
3115 
3116  end if
3117 
3118  !---------------------------------------------------------
3119 
3120  do it = 1, ntying
3121 
3122  do ip = 1, npoints_tying(it)
3123 
3124  !-------------------------------------------------
3125 
3126  naturalcoord(1) = tpcoord(ip, 1, it)
3127  naturalcoord(2) = tpcoord(ip, 2, it)
3128 
3129  call getshapefunc(fetype, naturalcoord, shapefunc)
3130 
3131  call getshapederiv(fetype, naturalcoord, shapederiv)
3132 
3133  !-------------------------------------------------
3134 
3135  do na = 1, nn
3136 
3137  do i = 1, 3
3138 
3139  u_rot(i, na) &
3140  = shapefunc(na) &
3141  *( zeta_ly*a_over_2_v3(i, na) )
3142 
3143  dudxi_rot(i, na) &
3144  = shapederiv(na, 1) &
3145  *( zeta_ly*a_over_2_v3(i, na) )
3146  dudeta_rot(i, na) &
3147  = shapederiv(na, 2) &
3148  *( zeta_ly*a_over_2_v3(i, na) )
3149  dudzeta_rot(i, na) &
3150  = shapefunc(na) &
3151  *( a_over_2_v3(i, na) )
3152 
3153  end do
3154 
3155  end do
3156 
3157  !-------------------------------------------------
3158 
3159  ! Covariant basis vector
3160  do i = 1, 3
3161 
3162  g1(i) = 0.0d0
3163  g2(i) = 0.0d0
3164  g3(i) = 0.0d0
3165 
3166  do na = 1, nn
3167 
3168  g1(i) = g1(i)+shapederiv(na, 1) &
3169  *elem(i, na) &
3170  +dudxi_rot(i, na)
3171  g2(i) = g2(i)+shapederiv(na, 2) &
3172  *elem(i, na) &
3173  +dudeta_rot(i, na)
3174  g3(i) = g3(i)+dudzeta_rot(i, na)
3175 
3176  end do
3177 
3178  end do
3179 
3180  !---------------------------------------------
3181 
3182  do i = 1, 3
3183 
3184  dudxi(i) = 0.0d0
3185  dudeta(i) = 0.0d0
3186  dudzeta(i) = 0.0d0
3187  dudxi_trans(i) = 0.0d0
3188  dudeta_trans(i) = 0.0d0
3189 
3190  do na = 1, nn
3191 
3192  dudxi(i) &
3193  = dudxi(i) &
3194  +shapederiv(na, 1) &
3195  *( edisp(i, na) &
3196  +zeta_ly*a_over_2_theta_cross_v3(i, na) )
3197  dudeta(i) &
3198  = dudeta(i) &
3199  +shapederiv(na, 2) &
3200  *( edisp(i, na) &
3201  +zeta_ly*a_over_2_theta_cross_v3(i, na) )
3202  dudzeta(i) &
3203  = dudzeta(i) &
3204  +shapefunc(na) &
3205  *( a_over_2_theta_cross_v3(i, na) )
3206  dudxi_trans(i) = dudxi_trans(i)+shapederiv(na, 1)*edisp(i, na)
3207  dudeta_trans(i) = dudeta_trans(i)+shapederiv(na, 2)*edisp(i, na)
3208 
3209  end do
3210 
3211  end do
3212 
3213  !---------------------------------------------
3214 
3215  if( use_gl_strain ) then
3216  g1_cur(:) = g1(:)+dudxi_trans(:)
3217  g2_cur(:) = g2(:)+dudeta_trans(:)
3218  g3_cur(:) = g3(:)
3219  g1_ref(:) = g1_cur(:)-dudxi(:)
3220  g2_ref(:) = g2_cur(:)-dudeta(:)
3221  g3_ref(:) = g3_cur(:)-dudzeta(:)
3222 
3223  e11_di(ip, it) = 0.5d0*(dot_product(g1_cur, g1_cur)-dot_product(g1_ref, g1_ref))
3224  e22_di(ip, it) = 0.5d0*(dot_product(g2_cur, g2_cur)-dot_product(g2_ref, g2_ref))
3225  e12_2 = dot_product(g1_cur, g2_cur)-dot_product(g1_ref, g2_ref)
3226  e23_di_2(ip, it) = dot_product(g2_cur, g3_cur)-dot_product(g2_ref, g3_ref)
3227  e31_di_2(ip, it) = dot_product(g3_cur, g1_cur)-dot_product(g3_ref, g1_ref)
3228  else
3229  ! Infinitesimal strain tensor
3230  e11_di(ip, it) &
3231  = 0.5d0 &
3232  *( ( g1(1)*dudxi(1) +dudxi(1) *g1(1) ) &
3233  +( g1(2)*dudxi(2) +dudxi(2) *g1(2) ) &
3234  +( g1(3)*dudxi(3) +dudxi(3) *g1(3) ) )
3235  e22_di(ip, it) &
3236  = 0.5d0 &
3237  *( ( g2(1)*dudeta(1)+dudeta(1)*g2(1) ) &
3238  +( g2(2)*dudeta(2)+dudeta(2)*g2(2) ) &
3239  +( g2(3)*dudeta(3)+dudeta(3)*g2(3) ) )
3240  e12_2 &
3241  = ( g1(1)*dudeta(1) +dudxi(1) *g2(1) ) &
3242  +( g1(2)*dudeta(2) +dudxi(2) *g2(2) ) &
3243  +( g1(3)*dudeta(3) +dudxi(3) *g2(3) )
3244  e23_di_2(ip, it) &
3245  = ( g2(1)*dudzeta(1)+dudeta(1) *g3(1) ) &
3246  +( g2(2)*dudzeta(2)+dudeta(2) *g3(2) ) &
3247  +( g2(3)*dudzeta(3)+dudeta(3) *g3(3) )
3248  e31_di_2(ip, it) &
3249  = ( g3(1)*dudxi(1) +dudzeta(1)*g1(1) ) &
3250  +( g3(2)*dudxi(2) +dudzeta(2)*g1(2) ) &
3251  +( g3(3)*dudxi(3) +dudzeta(3)*g1(3) )
3252  endif
3253 
3254  !-------------------------------------------------
3255 
3256  end do
3257 
3258  end do
3259 
3260  !--------------------------------------------------------
3261 
3262  call fstr_shell_layer_zeta( gausses(1), n_layer, zeta, zeta_ly, ierr_quad )
3263  if( ierr_quad /= 0 ) stop "Invalid shell layer zeta"
3264 
3265  !--------------------------------------------------------
3266 
3267  npoints = nn
3268  if( use_surface_gauss ) npoints = numofquadpoints(fetype)
3269 
3270  do lx = 1, npoints
3271 
3272  !--------------------------------------------------
3273 
3274  if( use_surface_gauss ) then
3275  call getquadpoint(fetype, lx, naturalcoord)
3276  else
3277  naturalcoord(1) = nncoord(lx, 1)
3278  naturalcoord(2) = nncoord(lx, 2)
3279  endif
3280 
3281  xi_lx = naturalcoord(1)
3282  eta_lx = naturalcoord(2)
3283 
3284  call getshapefunc(fetype, naturalcoord, shapefunc)
3285 
3286  call getshapederiv(fetype, naturalcoord, shapederiv)
3287 
3288  !--------------------------------------------------
3289 
3290  do na = 1, nn
3291 
3292  do i = 1, 3
3293 
3294  u_rot(i, na) &
3295  = shapefunc(na) &
3296  *( zeta_ly*a_over_2_v3(i, na) )
3297 
3298  dudxi_rot(i, na) &
3299  = shapederiv(na, 1) &
3300  *( zeta_ly*a_over_2_v3(i, na) )
3301  dudeta_rot(i, na) &
3302  = shapederiv(na, 2) &
3303  *( zeta_ly*a_over_2_v3(i, na) )
3304  dudzeta_rot(i, na) &
3305  = shapefunc(na) &
3306  *( a_over_2_v3(i, na) )
3307 
3308  end do
3309 
3310  end do
3311 
3312  !--------------------------------------------------
3313 
3314  ! Covariant basis vector
3315  do i = 1, 3
3316 
3317  g1(i) = 0.0d0
3318  g2(i) = 0.0d0
3319  g3(i) = 0.0d0
3320 
3321  do na = 1, nn
3322 
3323  g1(i) = g1(i)+shapederiv(na, 1) &
3324  *elem(i, na) &
3325  +dudxi_rot(i, na)
3326  g2(i) = g2(i)+shapederiv(na, 2) &
3327  *elem(i, na) &
3328  +dudeta_rot(i, na)
3329  g3(i) = g3(i)+dudzeta_rot(i, na)
3330 
3331  end do
3332 
3333  end do
3334 
3335  !--------------------------------------------------
3336 
3337  if( use_gl_strain ) then
3338  dudxi_trans(:) = 0.0d0
3339  dudeta_trans(:) = 0.0d0
3340  do na = 1, nn
3341  dudxi_trans(:) = dudxi_trans(:)+shapederiv(na, 1)*edisp(1:3, na)
3342  dudeta_trans(:) = dudeta_trans(:)+shapederiv(na, 2)*edisp(1:3, na)
3343  end do
3344  g1(:) = g1(:)+dudxi_trans(:)
3345  g2(:) = g2(:)+dudeta_trans(:)
3346  endif
3347 
3348  !--------------------------------------------------
3349 
3350  ! Jacobian
3351  det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
3352  +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
3353  +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )
3354  det_cur = det
3355 
3356  if(det == 0.0d0) then
3357  write(*,*)"ERROR:LIB Shell in l2009 Not Jacobian"
3358  stop
3359  endif
3360 
3361  det_inv = 1.0d0/det
3362 
3363  !--------------------------------------------------
3364 
3365  ! Contravariant basis vector
3366  cg1(1) = det_inv &
3367  *( g2(2)*g3(3)-g2(3)*g3(2) )
3368  cg1(2) = det_inv &
3369  *( g2(3)*g3(1)-g2(1)*g3(3) )
3370  cg1(3) = det_inv &
3371  *( g2(1)*g3(2)-g2(2)*g3(1) )
3372  cg2(1) = det_inv &
3373  *( g3(2)*g1(3)-g3(3)*g1(2) )
3374  cg2(2) = det_inv &
3375  *( g3(3)*g1(1)-g3(1)*g1(3) )
3376  cg2(3) = det_inv &
3377  *( g3(1)*g1(2)-g3(2)*g1(1) )
3378  cg3(1) = det_inv &
3379  *( g1(2)*g2(3)-g1(3)*g2(2) )
3380  cg3(2) = det_inv &
3381  *( g1(3)*g2(1)-g1(1)*g2(3) )
3382  cg3(3) = det_inv &
3383  *( g1(1)*g2(2)-g1(2)*g2(1) )
3384 
3385  !--------------------------------------------------
3386 
3387  g3_abs = dsqrt( g3(1)*g3(1) &
3388  +g3(2)*g3(2) &
3389  +g3(3)*g3(3) )
3390 
3391  !--------------------------------------------------
3392 
3393  ! Orthonormal vectors
3394 
3395  e3_hat(1) = g3(1)/g3_abs
3396  e3_hat(2) = g3(2)/g3_abs
3397  e3_hat(3) = g3(3)/g3_abs
3398 
3399  e1_hat(1) = g2(2)*e3_hat(3) &
3400  -g2(3)*e3_hat(2)
3401  e1_hat(2) = g2(3)*e3_hat(1) &
3402  -g2(1)*e3_hat(3)
3403  e1_hat(3) = g2(1)*e3_hat(2) &
3404  -g2(2)*e3_hat(1)
3405  e1_hat_abs = dsqrt( e1_hat(1)*e1_hat(1) &
3406  +e1_hat(2)*e1_hat(2) &
3407  +e1_hat(3)*e1_hat(3) )
3408  e1_hat(1) = e1_hat(1)/e1_hat_abs
3409  e1_hat(2) = e1_hat(2)/e1_hat_abs
3410  e1_hat(3) = e1_hat(3)/e1_hat_abs
3411 
3412  e2_hat(1) = e3_hat(2)*e1_hat(3) &
3413  -e3_hat(3)*e1_hat(2)
3414  e2_hat(2) = e3_hat(3)*e1_hat(1) &
3415  -e3_hat(1)*e1_hat(3)
3416  e2_hat(3) = e3_hat(1)*e1_hat(2) &
3417  -e3_hat(2)*e1_hat(1)
3418  e2_hat_abs = dsqrt( e2_hat(1)*e2_hat(1) &
3419  +e2_hat(2)*e2_hat(2) &
3420  +e2_hat(3)*e2_hat(3) )
3421  e2_hat(1) = e2_hat(1)/e2_hat_abs
3422  e2_hat(2) = e2_hat(2)/e2_hat_abs
3423  e2_hat(3) = e2_hat(3)/e2_hat_abs
3424 
3425  !--------------------------------------------------
3426 
3427  do i = 1, 3
3428 
3429  dudxi(i) = 0.0d0
3430  dudeta(i) = 0.0d0
3431  dudzeta(i) = 0.0d0
3432  dudxi_trans(i) = 0.0d0
3433  dudeta_trans(i) = 0.0d0
3434 
3435  do na = 1, nn
3436 
3437  dudxi(i) &
3438  = dudxi(i) &
3439  +shapederiv(na, 1) &
3440  *( edisp(i, na) &
3441  +zeta_ly*a_over_2_theta_cross_v3(i, na) )
3442  dudeta(i) &
3443  = dudeta(i) &
3444  +shapederiv(na, 2) &
3445  *( edisp(i, na) &
3446  +zeta_ly*a_over_2_theta_cross_v3(i, na) )
3447  dudzeta(i) &
3448  = dudzeta(i) &
3449  +shapefunc(na) &
3450  *( a_over_2_theta_cross_v3(i, na) )
3451  dudxi_trans(i) = dudxi_trans(i)+shapederiv(na, 1)*edisp(i, na)
3452  dudeta_trans(i) = dudeta_trans(i)+shapederiv(na, 2)*edisp(i, na)
3453 
3454  end do
3455 
3456  end do
3457 
3458  !--------------------------------------------------
3459 
3460  !--------------------------------------------------
3461 
3462  if( use_gl_strain ) then
3463  if( etype == fe_mitc4_shell ) then
3464  g1_cur(:) = g1(:)
3465  g2_cur(:) = g2(:)
3466  else
3467  g1_cur(:) = g1(:)+dudxi_trans(:)
3468  g2_cur(:) = g2(:)+dudeta_trans(:)
3469  endif
3470  g3_cur(:) = g3(:)
3471  g1_ref(:) = g1_cur(:)-dudxi(:)
3472  g2_ref(:) = g2_cur(:)-dudeta(:)
3473  g3_ref(:) = g3_cur(:)-dudzeta(:)
3474  det_ref = g1_ref(1)*( g2_ref(2)*g3_ref(3)-g2_ref(3)*g3_ref(2) ) &
3475  +g1_ref(2)*( g2_ref(3)*g3_ref(1)-g2_ref(1)*g3_ref(3) ) &
3476  +g1_ref(3)*( g2_ref(1)*g3_ref(2)-g2_ref(2)*g3_ref(1) )
3477 
3478  e11 = 0.5d0*(dot_product(g1_cur, g1_cur)-dot_product(g1_ref, g1_ref))
3479  e22 = 0.5d0*(dot_product(g2_cur, g2_cur)-dot_product(g2_ref, g2_ref))
3480  e12_2 = dot_product(g1_cur, g2_cur)-dot_product(g1_ref, g2_ref)
3481  e23_2 = dot_product(g2_cur, g3_cur)-dot_product(g2_ref, g3_ref)
3482  e31_2 = dot_product(g3_cur, g1_cur)-dot_product(g3_ref, g1_ref)
3483  else
3484  ! Infinitesimal strain tensor
3485  e11 &
3486  = 0.5d0 &
3487  *( ( g1(1)*dudxi(1) +dudxi(1) *g1(1) ) &
3488  +( g1(2)*dudxi(2) +dudxi(2) *g1(2) ) &
3489  +( g1(3)*dudxi(3) +dudxi(3) *g1(3) ) )
3490  e22 &
3491  = 0.5d0 &
3492  *( ( g2(1)*dudeta(1)+dudeta(1)*g2(1) ) &
3493  +( g2(2)*dudeta(2)+dudeta(2)*g2(2) ) &
3494  +( g2(3)*dudeta(3)+dudeta(3)*g2(3) ) )
3495  e12_2 &
3496  = ( g1(1)*dudeta(1) +dudxi(1) *g2(1) ) &
3497  +( g1(2)*dudeta(2) +dudxi(2) *g2(2) ) &
3498  +( g1(3)*dudeta(3) +dudxi(3) *g2(3) )
3499  e23_2 &
3500  = ( g2(1)*dudzeta(1)+dudeta(1) *g3(1) ) &
3501  +( g2(2)*dudzeta(2)+dudeta(2) *g3(2) ) &
3502  +( g2(3)*dudzeta(3)+dudeta(3) *g3(3) )
3503  e31_2 &
3504  = ( g3(1)*dudxi(1) +dudzeta(1)*g1(1) ) &
3505  +( g3(2)*dudxi(2) +dudzeta(2)*g1(2) ) &
3506  +( g3(3)*dudxi(3) +dudzeta(3)*g1(3) )
3507  endif
3508 
3509  ! MITC4
3510  if( etype .EQ. fe_mitc4_shell ) then
3511 
3512  e23_2 = 0.0d0
3513  e31_2 = 0.0d0
3514 
3515  ! e23_as_2
3516  e23_2 &
3517  = 0.5d0*( 1.0d0-xi_lx )*e23_di_2(4, 1) &
3518  +0.5d0*( 1.0d0+xi_lx )*e23_di_2(2, 1)
3519  ! e31_as_2
3520  e31_2 &
3521  = 0.5d0*( 1.0d0-eta_lx )*e31_di_2(1, 1) &
3522  +0.5d0*( 1.0d0+eta_lx )*e31_di_2(3, 1)
3523 
3524  ! MITC9
3525  else if( etype .EQ. fe_mitc9_shell ) then
3526 
3527  xxi_lx = xi_lx /dsqrt( 1.0d0/3.0d0 )
3528  eeta_lx = eta_lx/dsqrt( 3.0d0/5.0d0 )
3529 
3530  do ip = 1, npoints_tying(1)
3531 
3532  h(ip, 1) &
3533  = ( 0.5d0*( 1.0d0+xxi_di(ip, 1)*xxi_lx ) ) &
3534  *( ( 0.5d0*eeta_di(ip, 1)*eeta_lx ) &
3535  *( 1.0d0+eeta_di(ip, 1)*eeta_lx ) &
3536  +( 1.0d0-eeta_di(ip, 1)*eeta_di(ip, 1) ) &
3537  *( 1.0d0-eeta_lx*eeta_lx ) )
3538 
3539  end do
3540 
3541  xxi_lx = xi_lx /dsqrt( 3.0d0/5.0d0 )
3542  eeta_lx = eta_lx/dsqrt( 1.0d0/3.0d0 )
3543 
3544  do ip = 1, npoints_tying(2)
3545 
3546  h(ip, 2) &
3547  = ( ( 0.5d0*xxi_di(ip, 2) *xxi_lx ) &
3548  *( 1.0d0+xxi_di(ip, 2) *xxi_lx ) &
3549  +( 1.0d0-xxi_di(ip, 2) *xxi_di(ip, 2) ) &
3550  *( 1.0d0-xxi_lx*xxi_lx ) ) &
3551  *( 0.5d0*( 1.0d0+eeta_di(ip, 2)*eeta_lx ) )
3552 
3553  end do
3554 
3555  xxi_lx = xi_lx /dsqrt( 1.0d0/3.0d0 )
3556  eeta_lx = eta_lx/dsqrt( 1.0d0/3.0d0 )
3557 
3558  do ip = 1, npoints_tying(3)
3559 
3560  h(ip, 3) &
3561  = ( 0.5d0*( 1.0d0+xxi_di(ip, 1)*xxi_lx ) ) &
3562  *( 0.5d0*( 1.0d0+eeta_di(ip, 1)*eeta_lx ) )
3563 
3564  end do
3565 
3566  e11 = 0.0d0
3567  e31_2 = 0.0d0
3568 
3569  ! e11_as, e31_as_2
3570  do ip = 1, npoints_tying(1)
3571 
3572  e11 = e11 +h(ip, 1)*e11_di(ip, 1)
3573  e31_2 = e31_2+h(ip, 1)*e31_di_2(ip, 1)
3574 
3575  end do
3576 
3577  e22 = 0.0d0
3578  e23_2 = 0.0d0
3579 
3580  ! e22_as, e23_as_2
3581  do ip = 1, npoints_tying(2)
3582 
3583  e22 = e22 +h(ip, 2)*e22_di(ip, 2)
3584  e23_2 = e23_2+h(ip, 2)*e23_di_2(ip, 2)
3585 
3586  end do
3587 
3588  e12_2 = 0.0d0
3589 
3590  ! e12_as_2
3591  do ip = 1, npoints_tying(3)
3592 
3593  e12_2 = e12_2+h(ip, 3)*e12_di_2(ip, 3)
3594 
3595  end do
3596 
3597  ! MITC3
3598  else if( etype .EQ. fe_mitc3_shell ) then
3599 
3600  e23_2 = 0.0d0
3601  e31_2 = 0.0d0
3602 
3603  ! e23_as_2
3604  e23_2 &
3605  = ( 1.0d0-xi_lx )*e23_di_2(2, 1) &
3606  +xi_lx *e31_di_2(1, 1) &
3607  +xi_lx *( e23_di_2(3, 1)-e31_di_2(3, 1) )
3608  ! e31_as_2
3609  e31_2 &
3610  = eta_lx*e23_di_2(2, 1) &
3611  +( 1.0d0-eta_lx )*e31_di_2(1, 1) &
3612  -eta_lx*( e23_di_2(3, 1)-e31_di_2(3, 1) )
3613 
3614  end if
3615 
3616  !--------------------------------------------------
3617 
3618  ! { E } vector
3619  ev(1) = e11
3620  ev(2) = e22
3621  ev(3) = e12_2
3622  ev(4) = e23_2
3623  ev(5) = e31_2
3624 
3625  ! Infinitesimal strain tensor
3626  ! [ E ] matrix
3627  e(1, 1) = ev(1)
3628  e(2, 2) = ev(2)
3629  e(3, 3) = 0.0d0
3630  e(1, 2) = 0.5d0*ev(3)
3631  e(2, 1) = 0.5d0*ev(3)
3632  e(2, 3) = 0.5d0*ev(4)
3633  e(3, 2) = 0.5d0*ev(4)
3634  e(3, 1) = 0.5d0*ev(5)
3635  e(1, 3) = 0.5d0*ev(5)
3636 
3637  !--------------------------------------------------
3638  ! write(*,*) 'Stress_n_layer', n_layer
3639  e1_hat_mat(:) = e1_hat(:)
3640  e2_hat_mat(:) = e2_hat(:)
3641  e3_hat_mat(:) = e3_hat(:)
3642  cg1_mat(:) = cg1(:)
3643  cg2_mat(:) = cg2(:)
3644  cg3_mat(:) = cg3(:)
3645  if( use_gl_strain ) then
3646  call shell_basis_from_covariant(g1_ref, g2_ref, g3_ref, &
3647  e1_hat_mat, e2_hat_mat, e3_hat_mat, cg1_mat, cg2_mat, cg3_mat, det_mat)
3648  endif
3649 
3650  call matlmatrix_shell &
3651  (gausses(lx), shell, d, &
3652  e1_hat_mat, e2_hat_mat, e3_hat_mat, cg1_mat, cg2_mat, cg3_mat, &
3653  alpha, n_layer)
3654 
3655  !--------------------------------------------------
3656 
3657  sv = matmul( d, ev )
3658 
3659  if( present( local_stress_override ) ) then
3660  if( lx <= size(local_stress_override, 1) .and. size(local_stress_override, 2) >= 6 ) then
3661  sv(1) = local_stress_override(lx, 1)
3662  sv(2) = local_stress_override(lx, 2)
3663  sv(3) = local_stress_override(lx, 4)
3664  sv(4) = local_stress_override(lx, 5)
3665  sv(5) = local_stress_override(lx, 6)
3666  endif
3667  endif
3668 
3669  ! Infinitesimal stress tensor
3670  ! [ S ] matrix
3671  s(1, 1) = sv(1)
3672  s(2, 2) = sv(2)
3673  s(3, 3) = 0.0d0
3674  s(1, 2) = sv(3)
3675  s(2, 1) = sv(3)
3676  s(2, 3) = sv(4)
3677  s(3, 2) = sv(4)
3678  s(3, 1) = sv(5)
3679  s(1, 3) = sv(5)
3680 
3681  !------------------------------------------------
3682 
3683  if( present( local_strain ) ) then
3684  local_strain(lx, 1) = ev(1)
3685  local_strain(lx, 2) = ev(2)
3686  local_strain(lx, 3) = 0.0d0
3687  local_strain(lx, 4) = ev(3)
3688  local_strain(lx, 5) = ev(4)
3689  local_strain(lx, 6) = ev(5)
3690  endif
3691 
3692  if( present( local_stress ) ) then
3693  local_stress(lx, 1) = sv(1)
3694  local_stress(lx, 2) = sv(2)
3695  local_stress(lx, 3) = 0.0d0
3696  local_stress(lx, 4) = sv(3)
3697  local_stress(lx, 5) = sv(4)
3698  local_stress(lx, 6) = sv(5)
3699  endif
3700 
3701  if( use_gl_strain ) then
3702  g1(:) = g1_ref(:)
3703  g2(:) = g2_ref(:)
3704  g3(:) = g3_ref(:)
3705 
3706  det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
3707  +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
3708  +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )
3709  if(det == 0.0d0) then
3710  write(*,*)"ERROR:LIB Shell in l2009 Not Jacobian"
3711  stop
3712  endif
3713 
3714  det_inv = 1.0d0/det
3715  cg1(1) = det_inv*( g2(2)*g3(3)-g2(3)*g3(2) )
3716  cg1(2) = det_inv*( g2(3)*g3(1)-g2(1)*g3(3) )
3717  cg1(3) = det_inv*( g2(1)*g3(2)-g2(2)*g3(1) )
3718  cg2(1) = det_inv*( g3(2)*g1(3)-g3(3)*g1(2) )
3719  cg2(2) = det_inv*( g3(3)*g1(1)-g3(1)*g1(3) )
3720  cg2(3) = det_inv*( g3(1)*g1(2)-g3(2)*g1(1) )
3721  cg3(1) = det_inv*( g1(2)*g2(3)-g1(3)*g2(2) )
3722  cg3(2) = det_inv*( g1(3)*g2(1)-g1(1)*g2(3) )
3723  cg3(3) = det_inv*( g1(1)*g2(2)-g1(2)*g2(1) )
3724  cg1_ref(:) = cg1(:)
3725  cg2_ref(:) = cg2(:)
3726  cg3_ref(:) = cg3(:)
3727  endif
3728 
3729  stress(lx,1) &
3730  = ( s(1, 1)*g1(1)*g1(1) &
3731  +s(1, 2)*g1(1)*g2(1) &
3732  +s(1, 3)*g1(1)*g3(1) &
3733  +s(2, 1)*g2(1)*g1(1) &
3734  +s(2, 2)*g2(1)*g2(1) &
3735  +s(2, 3)*g2(1)*g3(1) &
3736  +s(3, 1)*g3(1)*g1(1) &
3737  +s(3, 2)*g3(1)*g2(1) )
3738  stress(lx,2) &
3739  = ( s(1, 1)*g1(2)*g1(2) &
3740  +s(1, 2)*g1(2)*g2(2) &
3741  +s(1, 3)*g1(2)*g3(2) &
3742  +s(2, 1)*g2(2)*g1(2) &
3743  +s(2, 2)*g2(2)*g2(2) &
3744  +s(2, 3)*g2(2)*g3(2) &
3745  +s(3, 1)*g3(2)*g1(2) &
3746  +s(3, 2)*g3(2)*g2(2) )
3747  stress(lx,3) &
3748  = ( s(1, 1)*g1(3)*g1(3) &
3749  +s(1, 2)*g1(3)*g2(3) &
3750  +s(1, 3)*g1(3)*g3(3) &
3751  +s(2, 1)*g2(3)*g1(3) &
3752  +s(2, 2)*g2(3)*g2(3) &
3753  +s(2, 3)*g2(3)*g3(3) &
3754  +s(3, 1)*g3(3)*g1(3) &
3755  +s(3, 2)*g3(3)*g2(3) )
3756  stress(lx,4) &
3757  = ( s(1, 1)*g1(1)*g1(2) &
3758  +s(1, 2)*g1(1)*g2(2) &
3759  +s(1, 3)*g1(1)*g3(2) &
3760  +s(2, 1)*g2(1)*g1(2) &
3761  +s(2, 2)*g2(1)*g2(2) &
3762  +s(2, 3)*g2(1)*g3(2) &
3763  +s(3, 1)*g3(1)*g1(2) &
3764  +s(3, 2)*g3(1)*g2(2) )
3765  stress(lx,5) &
3766  = ( s(1, 1)*g1(2)*g1(3) &
3767  +s(1, 2)*g1(2)*g2(3) &
3768  +s(1, 3)*g1(2)*g3(3) &
3769  +s(2, 1)*g2(2)*g1(3) &
3770  +s(2, 2)*g2(2)*g2(3) &
3771  +s(2, 3)*g2(2)*g3(3) &
3772  +s(3, 1)*g3(2)*g1(3) &
3773  +s(3, 2)*g3(2)*g2(3) )
3774  stress(lx,6) &
3775  = ( s(1, 1)*g1(3)*g1(1) &
3776  +s(1, 2)*g1(3)*g2(1) &
3777  +s(1, 3)*g1(3)*g3(1) &
3778  +s(2, 1)*g2(3)*g1(1) &
3779  +s(2, 2)*g2(3)*g2(1) &
3780  +s(2, 3)*g2(3)*g3(1) &
3781  +s(3, 1)*g3(3)*g1(1) &
3782  +s(3, 2)*g3(3)*g2(1) )
3783 
3784  strain(lx,1) &
3785  = ( e(1, 1)*cg1(1)*cg1(1) &
3786  +e(1, 2)*cg1(1)*cg2(1) &
3787  +e(1, 3)*cg1(1)*cg3(1) &
3788  +e(2, 1)*cg2(1)*cg1(1) &
3789  +e(2, 2)*cg2(1)*cg2(1) &
3790  +e(2, 3)*cg2(1)*cg3(1) &
3791  +e(3, 1)*cg3(1)*cg1(1) &
3792  +e(3, 2)*cg3(1)*cg2(1) )
3793  strain(lx,2) &
3794  = ( e(1, 1)*cg1(2)*cg1(2) &
3795  +e(1, 2)*cg1(2)*cg2(2) &
3796  +e(1, 3)*cg1(2)*cg3(2) &
3797  +e(2, 1)*cg2(2)*cg1(2) &
3798  +e(2, 2)*cg2(2)*cg2(2) &
3799  +e(2, 3)*cg2(2)*cg3(2) &
3800  +e(3, 1)*cg3(2)*cg1(2) &
3801  +e(3, 2)*cg3(2)*cg2(2) )
3802  strain(lx,3) &
3803  = ( e(1, 1)*cg1(3)*cg1(3) &
3804  +e(1, 2)*cg1(3)*cg2(3) &
3805  +e(1, 3)*cg1(3)*cg3(3) &
3806  +e(2, 1)*cg2(3)*cg1(3) &
3807  +e(2, 2)*cg2(3)*cg2(3) &
3808  +e(2, 3)*cg2(3)*cg3(3) &
3809  +e(3, 1)*cg3(3)*cg1(3) &
3810  +e(3, 2)*cg3(3)*cg2(3) )
3811  strain(lx,4) &
3812  = ( e(1, 1)*cg1(1)*cg1(2) &
3813  +e(1, 2)*cg1(1)*cg2(2) &
3814  +e(1, 3)*cg1(1)*cg3(2) &
3815  +e(2, 1)*cg2(1)*cg1(2) &
3816  +e(2, 2)*cg2(1)*cg2(2) &
3817  +e(2, 3)*cg2(1)*cg3(2) &
3818  +e(3, 1)*cg3(1)*cg1(2) &
3819  +e(3, 2)*cg3(1)*cg2(2) )
3820  strain(lx,5) &
3821  = ( e(1, 1)*cg1(2)*cg1(3) &
3822  +e(1, 2)*cg1(2)*cg2(3) &
3823  +e(1, 3)*cg1(2)*cg3(3) &
3824  +e(2, 1)*cg2(2)*cg1(3) &
3825  +e(2, 2)*cg2(2)*cg2(3) &
3826  +e(2, 3)*cg2(2)*cg3(3) &
3827  +e(3, 1)*cg3(2)*cg1(3) &
3828  +e(3, 2)*cg3(2)*cg2(3) )
3829  strain(lx,6) &
3830  = ( e(1, 1)*cg1(3)*cg1(1) &
3831  +e(1, 2)*cg1(3)*cg2(1) &
3832  +e(1, 3)*cg1(3)*cg3(1) &
3833  +e(2, 1)*cg2(3)*cg1(1) &
3834  +e(2, 2)*cg2(3)*cg2(1) &
3835  +e(2, 3)*cg2(3)*cg3(1) &
3836  +e(3, 1)*cg3(3)*cg1(1) &
3837  +e(3, 2)*cg3(3)*cg2(1) )
3838 
3839  if( use_gl_strain .and. opsstype == kopss_solution ) then
3840  jac = det_cur/det_ref
3841  if( jac == 0.0d0 ) stop "Fail to convert shell stress: detF=0"
3842 
3843  cauchy(:, :) = 0.0d0
3844  cauchy(:, :) = cauchy(:, :) &
3845  +s(1, 1)*outer_product3(g1_cur, g1_cur) &
3846  +s(1, 2)*outer_product3(g1_cur, g2_cur) &
3847  +s(1, 3)*outer_product3(g1_cur, g3_cur) &
3848  +s(2, 1)*outer_product3(g2_cur, g1_cur) &
3849  +s(2, 2)*outer_product3(g2_cur, g2_cur) &
3850  +s(2, 3)*outer_product3(g2_cur, g3_cur) &
3851  +s(3, 1)*outer_product3(g3_cur, g1_cur) &
3852  +s(3, 2)*outer_product3(g3_cur, g2_cur) &
3853  +s(3, 3)*outer_product3(g3_cur, g3_cur)
3854  cauchy(:, :) = cauchy(:, :)/jac
3855 
3856  cg_metric(1, 1) = dot_product(cg1_ref, cg1_ref)
3857  cg_metric(1, 2) = dot_product(cg1_ref, cg2_ref)
3858  cg_metric(1, 3) = dot_product(cg1_ref, cg3_ref)
3859  cg_metric(2, 1) = cg_metric(1, 2)
3860  cg_metric(2, 2) = dot_product(cg2_ref, cg2_ref)
3861  cg_metric(2, 3) = dot_product(cg2_ref, cg3_ref)
3862  cg_metric(3, 1) = cg_metric(1, 3)
3863  cg_metric(3, 2) = cg_metric(2, 3)
3864  cg_metric(3, 3) = dot_product(cg3_ref, cg3_ref)
3865 
3866  stretch_b(:, :) = 0.0d0
3867  stretch_b(:, :) = stretch_b(:, :) &
3868  +cg_metric(1, 1)*outer_product3(g1_cur, g1_cur) &
3869  +cg_metric(1, 2)*outer_product3(g1_cur, g2_cur) &
3870  +cg_metric(1, 3)*outer_product3(g1_cur, g3_cur) &
3871  +cg_metric(2, 1)*outer_product3(g2_cur, g1_cur) &
3872  +cg_metric(2, 2)*outer_product3(g2_cur, g2_cur) &
3873  +cg_metric(2, 3)*outer_product3(g2_cur, g3_cur) &
3874  +cg_metric(3, 1)*outer_product3(g3_cur, g1_cur) &
3875  +cg_metric(3, 2)*outer_product3(g3_cur, g2_cur) &
3876  +cg_metric(3, 3)*outer_product3(g3_cur, g3_cur)
3877 
3878  tensor(1) = stretch_b(1, 1)
3879  tensor(2) = stretch_b(2, 2)
3880  tensor(3) = stretch_b(3, 3)
3881  tensor(4) = stretch_b(1, 2)
3882  tensor(5) = stretch_b(2, 3)
3883  tensor(6) = stretch_b(3, 1)
3884  call get_principal(tensor, eigval, princ)
3885 
3886  do i = 1, 3
3887  if( eigval(i) <= 0.0d0 ) stop "Fail to calc shell log strain: stretch<0"
3888  eigval(i) = 0.5d0*dlog(eigval(i))
3889  eig_norm = dsqrt(dot_product(princ(1:3, i), princ(1:3, i)))
3890  if( eig_norm <= 0.0d0 ) stop "Fail to calc shell log strain: direction vector=0"
3891  princ(1:3, i) = princ(1:3, i)/eig_norm
3892  end do
3893 
3894  logstrain(:, :) = 0.0d0
3895  do i = 1, 3
3896  logstrain(:, :) = logstrain(:, :) &
3897  +eigval(i)*outer_product3(princ(1:3, i), princ(1:3, i))
3898  end do
3899 
3900  stress(lx, 1) = cauchy(1, 1)
3901  stress(lx, 2) = cauchy(2, 2)
3902  stress(lx, 3) = cauchy(3, 3)
3903  stress(lx, 4) = cauchy(1, 2)
3904  stress(lx, 5) = cauchy(2, 3)
3905  stress(lx, 6) = cauchy(3, 1)
3906 
3907  strain(lx, 1) = logstrain(1, 1)
3908  strain(lx, 2) = logstrain(2, 2)
3909  strain(lx, 3) = logstrain(3, 3)
3910  strain(lx, 4) = 2.0d0*logstrain(1, 2)
3911  strain(lx, 5) = 2.0d0*logstrain(2, 3)
3912  strain(lx, 6) = 2.0d0*logstrain(3, 1)
3913  endif
3914 
3915  !--------------------------------------------------
3916 
3917  end do
3918  !--------------------------------------------------------------------
3919 
3920  return
3921 
3922  !####################################################################
3923  end subroutine elementstress_shell_mitc
3924  !####################################################################
3925 
3926  !####################################################################
3927  subroutine updatestressshellul_elastic( gauss, dstrain, dstress, trace_coeff )
3928  !####################################################################
3929 
3930  use mmechgauss
3931 
3932  !--------------------------------------------------------------------
3933 
3934  type(tgaussstatus), intent(inout) :: gauss
3935  real(kind = kreal), intent(in) :: dstrain(6)
3936  real(kind = kreal), intent(in) :: dstress(6)
3937  real(kind = kreal), intent(in), optional :: trace_coeff
3938 
3939  !--------------------------------------------------------------------
3940 
3941  real(kind = kreal) :: dstress_obj(6)
3942 
3943  !--------------------------------------------------------------------
3944 
3945  call shellobjectivestressincrement( gauss%stress_bak(1:6), dstrain(1:6), dstress_obj, trace_coeff )
3946 
3947  gauss%strain(1:6) = gauss%strain_bak(1:6) + dstrain(1:6)
3948  gauss%stress(1:6) = gauss%stress_bak(1:6) + dstress_obj(1:6) + dstress(1:6)
3949 
3950  gauss%strain_energy = gauss%strain_energy_bak &
3951  + dot_product( gauss%stress(1:6), dstrain(1:6) )
3952  gauss%strain_energy = gauss%strain_energy &
3953  - 0.5d0*dot_product( dstress(1:6), dstrain(1:6) )
3954  gauss%strain_out(1:6) = gauss%strain(1:6)
3955  gauss%stress_out(1:6) = gauss%stress(1:6)
3956 
3957  return
3958 
3959  !####################################################################
3960  end subroutine updatestressshellul_elastic
3961  !####################################################################
3962 
3963  !####################################################################
3964  subroutine updateshelllayergauss_shell_mitc &
3965  (etype, nn, ndof, ecoord, element, edisp, thick, nddirector, ndrefdirector, ndbase_disp)
3966  !####################################################################
3967 
3968  use mmechgauss
3969  use quadrature
3970  use mmaterial, only: updatelag, iselastic
3971 
3972  !--------------------------------------------------------------------
3973 
3974  integer(kind = kint), intent(in) :: etype
3975  integer(kind = kint), intent(in) :: nn
3976  integer(kind = kint), intent(in) :: ndof
3977  real(kind = kreal), intent(in) :: ecoord(3, nn)
3978  type(telement), intent(inout) :: element
3979  real(kind = kreal), intent(in) :: edisp(6, nn)
3980  real(kind = kreal), intent(in) :: thick
3981  real(kind = kreal), intent(in), optional :: nddirector(3, nn)
3982  real(kind = kreal), intent(in), optional :: ndrefdirector(3, nn)
3983  real(kind = kreal), intent(in), optional :: ndbase_disp(6, nn)
3984 
3985  !--------------------------------------------------------------------
3986 
3987  integer(kind = kint) :: ng, ig, ilayer, ithick, ishell, ierr
3988  integer(kind = kint) :: flag
3989  real(kind = kreal) :: zeta, weight
3990  real(kind = kreal) :: trace_coeff
3991  real(kind = kreal) :: gpstrain(9, 6), gpstress(9, 6)
3992  real(kind = kreal) :: local_gpstrain(9, 6), local_gpstress(9, 6)
3993 
3994  !--------------------------------------------------------------------
3995 
3996  if( .not. associated( element%gausses ) ) return
3997  if( .not. associated( element%shell_layer_gausses ) ) return
3998  if( element%shell_nlayer <= 0 .or. element%shell_nthick <= 0 ) return
3999  if( .not. ( shellsupportsfiniterotationkinematics( etype, nn ) &
4000  .and. iselastic(element%gausses(1)%pMaterial%mtype) ) ) return
4001 
4002  flag = element%gausses(1)%pMaterial%nlgeom_flag
4003  ng = numofquadpoints( etype )
4004  if( ng <= 0 .or. ng > 9 ) return
4005 
4006  do ilayer = 1, element%shell_nlayer
4007  do ithick = 1, element%shell_nthick
4008  call fstr_shell_thickness_quadrature( etype, ithick, zeta, weight, ierr )
4009  if( ierr /= 0 ) cycle
4010 
4011  gpstrain(:, :) = 0.0d0
4012  gpstress(:, :) = 0.0d0
4013  local_gpstrain(:, :) = 0.0d0
4014  local_gpstress(:, :) = 0.0d0
4015  call elementstress_shell_mitc( etype, nn, ndof, ecoord, element%gausses, edisp, &
4016  gpstrain(1:ng, 1:6), gpstress(1:ng, 1:6), thick, zeta, ilayer, element%shell_nlayer, &
4017  surface_gauss_points=.true., local_strain=local_gpstrain(1:ng, 1:6), &
4018  local_stress=local_gpstress(1:ng, 1:6), nddirector=nddirector, ndrefdirector=ndrefdirector, &
4019  ndbase_disp=ndbase_disp )
4020 
4021  do ig = 1, ng
4022  ishell = fstr_shell_layer_gauss_index( element, ig, ilayer, ithick )
4023  if( ishell <= 0 ) cycle
4024  if( flag == updatelag ) then
4025  trace_coeff = shellplanestresstracecoeff(element%shell_layer_gausses(ishell), ilayer)
4026  call updatestressshellul_elastic( element%shell_layer_gausses(ishell), &
4027  local_gpstrain(ig, 1:6), local_gpstress(ig, 1:6), trace_coeff=trace_coeff )
4028  else
4029  element%shell_layer_gausses(ishell)%strain(1:6) = local_gpstrain(ig, 1:6)
4030  element%shell_layer_gausses(ishell)%stress(1:6) = local_gpstress(ig, 1:6)
4031  element%shell_layer_gausses(ishell)%strain_energy = &
4032  0.5d0*dot_product( local_gpstress(ig, 1:6), local_gpstrain(ig, 1:6) )
4033  element%shell_layer_gausses(ishell)%strain_out(1:6) = gpstrain(ig, 1:6)
4034  element%shell_layer_gausses(ishell)%stress_out(1:6) = gpstress(ig, 1:6)
4035  endif
4036  enddo
4037  enddo
4038  enddo
4039 
4040  return
4041 
4042  !####################################################################
4043  end subroutine updateshelllayergauss_shell_mitc
4044  !####################################################################
4045 
4046  !####################################################################
4047  subroutine dl_shell &
4048  (etype, nn, ndof, xx, yy, zz, rho, thick, &
4049  ltype, params, vect, nsize, gausses)
4050  !####################################################################
4051 
4052  use hecmw
4053  use m_utilities
4054  use mmechgauss
4055  use quadrature
4056 
4057  type(tgaussstatus), intent(in) :: gausses(:)
4058  !--------------------------------------------------------------------
4059 
4060  integer(kind = kint), intent(in) :: etype
4061  integer(kind = kint), intent(in) :: nn
4062  integer(kind = kint), intent(in) :: ndof
4063  real(kind = kreal), intent(in) :: xx(*), yy(*), zz(*)
4064  real(kind = kreal), intent(in) :: rho
4065  real(kind = kreal), intent(in) :: thick
4066  real(kind = kreal), intent(in) :: params(*)
4067  real(kind = kreal), intent(out) :: vect(*)
4068  integer(kind = kint), intent(out) :: nsize
4069 
4070  !--------------------------------------------------------------------
4071 
4072  integer :: ivol, isurf
4073  integer :: lx, ly
4074  integer :: fetype
4075  integer :: ny
4076  integer :: i
4077  integer(kind=kint) :: ierr_quad
4078  integer :: na, nb
4079  integer :: isize
4080  integer :: jsize1, jsize2, jsize3, &
4081  jsize4, jsize5, jsize6
4082  integer :: ltype
4083  integer :: n_totlyr, n_layer
4084 
4085  real(kind = kreal) :: elem(3, nn)
4086  real(kind = kreal) :: val
4087  real(kind = kreal) :: ax, ay, az
4088  real(kind = kreal) :: rx, ry, rz
4089  real(kind = kreal) :: xi_lx, eta_lx, zeta_ly
4090  real(kind = kreal) :: w_w_lx, w_ly
4091  real(kind = kreal) :: naturalcoord(2)
4092  real(kind = kreal) :: nncoord(nn, 2)
4093  real(kind = kreal) :: shapefunc(nn)
4094  real(kind = kreal) :: shapederiv(nn, 2)
4095  real(kind = kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
4096  real(kind = kreal) :: v1_abs, v2_abs, v3_abs
4097  real(kind = kreal) :: a_over_2_v3(3, nn)
4098  real(kind = kreal) :: u_rot(3, nn)
4099  real(kind = kreal) :: dudxi_rot(3, nn), dudeta_rot(3, nn), &
4100  dudzeta_rot(3, nn)
4101  real(kind = kreal) :: g1(3), g2(3), g3(3)
4102  real(kind = kreal) :: g1_cross_g2(3)
4103  real(kind = kreal) :: e_0(3)
4104  real(kind = kreal) :: det
4105  real(kind = kreal) :: det_cg3(3)
4106  real(kind = kreal) :: det_cg3_abs
4107  real(kind = kreal) :: w_w_w_det
4108  real(kind = kreal) :: n(3, ndof*nn)
4109  real(kind = kreal) :: hx, hy, hz
4110  real(kind = kreal) :: phx, phy, phz
4111  real(kind = kreal) :: coefx, coefy, coefz
4112  real(kind = kreal) :: x, y, z
4113 
4114  ny = 0
4115 
4116  !--------------------------------------------------------------------
4117 
4118  ! BX LTYPE=1 :BODY FORCE IN X-DIRECTION
4119  ! BY LTYPE=2 :BODY FORCE IN Y-DIRECTION
4120  ! BZ LTYPE=3 :BODY FORCE IN Z-DIRECTION
4121  ! CRAV LTYPE=4 :GRAVITY FORCE
4122  ! CENT LTYPE=5 :CENTRIFUGAL LOAD
4123  ! P LTYPE=10 :TRACTION IN NORMAL-DIRECTION FOR SHELL SURFACE
4124 
4125  !--------------------------------------------------------------------
4126 
4127  nsize = ndof*nn
4128 
4129  !--------------------------------------------------------------------
4130 
4131  val = params(1)
4132  ax = params(2)
4133  ay = params(3)
4134  az = params(4)
4135  rx = params(5)
4136  ry = params(6)
4137  rz = params(7)
4138 
4139  !--------------------------------------------------------------------
4140 
4141  ! MITC4
4142  if( etype .EQ. fe_mitc4_shell ) then
4143 
4144  fetype = fe_mitc4_shell
4145 
4146  ny = 2
4147 
4148  ! MITC9
4149  else if( etype .EQ. fe_mitc9_shell ) then
4150 
4151  fetype = fe_mitc9_shell
4152 
4153  ny = 3
4154 
4155  ! MITC3
4156  else if( etype .EQ. fe_mitc3_shell ) then
4157 
4158  fetype = fe_mitc3_shell
4159 
4160  ny = 2
4161 
4162  end if
4163 
4164  !--------------------------------------------------------------------
4165 
4166  do na = 1, nn
4167 
4168  elem(1, na) = xx(na)
4169  elem(2, na) = yy(na)
4170  elem(3, na) = zz(na)
4171 
4172  end do
4173 
4174  !-------------------------------------------------------------------
4175 
4176  ! xi-coordinate at a node in a local element
4177  ! eta-coordinate at a node in a local element
4178  call getnodalnaturalcoord(fetype, nncoord)
4179 
4180  !--------------------------------------------------------------------
4181 
4182  ! Local load vector
4183  do isize = 1, ndof*nn
4184 
4185  vect(isize) = 0.0d0
4186 
4187  end do
4188 
4189  !--------------------------------------------------------------------
4190 
4191  ! xi-coordinate at the center point in a local element
4192  ! eta-coordinate at the center point in a local element
4193  naturalcoord(1) = 0.0d0
4194  naturalcoord(2) = 0.0d0
4195 
4196  call getshapederiv(fetype, naturalcoord, shapederiv)
4197 
4198  !--------------------------------------------------------------
4199 
4200  ! Covariant basis vector
4201  do i = 1, 3
4202 
4203  g1(i) = 0.0d0
4204 
4205  do na = 1, nn
4206 
4207  g1(i) = g1(i)+shapederiv(na, 1) &
4208  *elem(i, na)
4209 
4210  end do
4211 
4212  end do
4213 
4214  e_0(1) = g1(1)
4215  e_0(2) = g1(2)
4216  e_0(3) = g1(3)
4217 
4218  !--------------------------------------------------------------
4219 
4220  do nb = 1, nn
4221 
4222  !--------------------------------------------------------
4223 
4224  naturalcoord(1) = nncoord(nb, 1)
4225  naturalcoord(2) = nncoord(nb, 2)
4226 
4227  call getshapederiv(fetype, naturalcoord, shapederiv)
4228 
4229  !--------------------------------------------------------
4230 
4231  ! Covariant basis vector
4232  do i = 1, 3
4233 
4234  g1(i) = 0.0d0
4235  g2(i) = 0.0d0
4236 
4237  do na = 1, nn
4238 
4239  g1(i) = g1(i)+shapederiv(na, 1) &
4240  *elem(i, na)
4241  g2(i) = g2(i)+shapederiv(na, 2) &
4242  *elem(i, na)
4243 
4244  end do
4245 
4246  end do
4247 
4248  !--------------------------------------------------------
4249 
4250  det_cg3(1) = g1(2)*g2(3)-g1(3)*g2(2)
4251  det_cg3(2) = g1(3)*g2(1)-g1(1)*g2(3)
4252  det_cg3(3) = g1(1)*g2(2)-g1(2)*g2(1)
4253 
4254  det_cg3_abs = dsqrt( det_cg3(1)*det_cg3(1) &
4255  +det_cg3(2)*det_cg3(2) &
4256  +det_cg3(3)*det_cg3(3) )
4257 
4258  v3(1, nb) = det_cg3(1)/det_cg3_abs
4259  v3(2, nb) = det_cg3(2)/det_cg3_abs
4260  v3(3, nb) = det_cg3(3)/det_cg3_abs
4261 
4262  !--------------------------------------------------------
4263 
4264  v2(1, nb) = v3(2, nb)*e_0(3)-v3(3, nb)*e_0(2)
4265  v2(2, nb) = v3(3, nb)*e_0(1)-v3(1, nb)*e_0(3)
4266  v2(3, nb) = v3(1, nb)*e_0(2)-v3(2, nb)*e_0(1)
4267 
4268  v2_abs = dsqrt( v2(1, nb)*v2(1, nb) &
4269  +v2(2, nb)*v2(2, nb) &
4270  +v2(3, nb)*v2(3, nb) )
4271 
4272  if( v2_abs .GT. 1.0d-15 ) then
4273 
4274  v2(1, nb) = v2(1, nb)/v2_abs
4275  v2(2, nb) = v2(2, nb)/v2_abs
4276  v2(3, nb) = v2(3, nb)/v2_abs
4277 
4278  v1(1, nb) = v2(2, nb)*v3(3, nb) &
4279  -v2(3, nb)*v3(2, nb)
4280  v1(2, nb) = v2(3, nb)*v3(1, nb) &
4281  -v2(1, nb)*v3(3, nb)
4282  v1(3, nb) = v2(1, nb)*v3(2, nb) &
4283  -v2(2, nb)*v3(1, nb)
4284 
4285  v1_abs = dsqrt( v1(1, nb)*v1(1, nb) &
4286  +v1(2, nb)*v1(2, nb) &
4287  +v1(3, nb)*v1(3, nb) )
4288 
4289  v1(1, nb) = v1(1, nb)/v1_abs
4290  v1(2, nb) = v1(2, nb)/v1_abs
4291  v1(3, nb) = v1(3, nb)/v1_abs
4292 
4293  else
4294 
4295  v1(1, nb) = 0.0d0
4296  v1(2, nb) = 0.0d0
4297  v1(3, nb) = -1.0d0
4298 
4299  v2(1, nb) = 0.0d0
4300  v2(2, nb) = 1.0d0
4301  v2(3, nb) = 0.0d0
4302 
4303  end if
4304 
4305  !--------------------------------------------------------
4306 
4307  v3(1, nb) = v1(2, nb)*v2(3, nb) &
4308  -v1(3, nb)*v2(2, nb)
4309  v3(2, nb) = v1(3, nb)*v2(1, nb) &
4310  -v1(1, nb)*v2(3, nb)
4311  v3(3, nb) = v1(1, nb)*v2(2, nb) &
4312  -v1(2, nb)*v2(1, nb)
4313 
4314  v3_abs = dsqrt( v3(1, nb)*v3(1, nb) &
4315  +v3(2, nb)*v3(2, nb) &
4316  +v3(3, nb)*v3(3, nb) )
4317 
4318  v3(1, nb) = v3(1, nb)/v3_abs
4319  v3(2, nb) = v3(2, nb)/v3_abs
4320  v3(3, nb) = v3(3, nb)/v3_abs
4321 
4322  !--------------------------------------------------------
4323 
4324  a_over_2_v3(1, nb) = 0.5d0*thick*v3(1, nb)
4325  a_over_2_v3(2, nb) = 0.5d0*thick*v3(2, nb)
4326  a_over_2_v3(3, nb) = 0.5d0*thick*v3(3, nb)
4327 
4328  !--------------------------------------------------------
4329 
4330  end do
4331 
4332  !--------------------------------------------------------------------
4333 
4334  ! Selection of load type
4335 
4336  ivol = 0
4337  isurf = 0
4338 
4339  if( ltype .LT. 10 ) then
4340 
4341  ivol = 1
4342 
4343  else if( ltype .GE. 10 ) then
4344 
4345  isurf = 1
4346 
4347  end if
4348 
4349  !--------------------------------------------------------------------
4350 
4351  !** Surface load
4352  if( isurf .EQ. 1 ) then
4353 
4354  !--------------------------------------------------------
4355 
4356  do lx = 1, numofquadpoints(fetype)
4357 
4358  !--------------------------------------------------
4359 
4360  call getquadpoint(fetype, lx, naturalcoord)
4361 
4362  xi_lx = naturalcoord(1)
4363  eta_lx = naturalcoord(2)
4364 
4365  w_w_lx = getweight(fetype, lx)
4366 
4367  call getshapefunc(fetype, naturalcoord, shapefunc)
4368 
4369  call getshapederiv(fetype, naturalcoord, shapederiv)
4370 
4371  !--------------------------------------------------
4372 
4373  do na = 1, nn
4374 
4375  do i = 1, 3
4376 
4377  u_rot(i, na) &
4378  = shapefunc(na) &
4379  *( 0.0d0*a_over_2_v3(i, na) )
4380 
4381  dudxi_rot(i, na) &
4382  = shapederiv(na, 1) &
4383  *( 0.0d0*a_over_2_v3(i, na) )
4384  dudeta_rot(i, na) &
4385  = shapederiv(na, 2) &
4386  *( 0.0d0*a_over_2_v3(i, na) )
4387  dudzeta_rot(i, na) &
4388  = shapefunc(na) &
4389  *( a_over_2_v3(i, na) )
4390 
4391  end do
4392 
4393  end do
4394 
4395  !--------------------------------------------------
4396 
4397  ! Covariant basis vector
4398  do i = 1, 3
4399 
4400  g1(i) = 0.0d0
4401  g2(i) = 0.0d0
4402  !g3(i) = 0.0D0
4403 
4404  do na = 1, nn
4405 
4406  g1(i) = g1(i)+shapederiv(na, 1) &
4407  *elem(i, na) &
4408  +dudxi_rot(i, na)
4409  g2(i) = g2(i)+shapederiv(na, 2) &
4410  *elem(i, na) &
4411  +dudeta_rot(i, na)
4412  !g3(i) = g3(i)+dudzeta_rot(i, na)
4413 
4414  end do
4415 
4416  end do
4417 
4418  !--------------------------------------------------
4419 
4420  !g3_abs = DSQRT( g3(1)*g3(1) &
4421  ! +g3(2)*g3(2) &
4422  ! +g3(3)*g3(3) )
4423 
4424  !--------------------------------------------------
4425 
4426  !e3_hat(1) = g3(1)/g3_abs
4427  !e3_hat(2) = g3(2)/g3_abs
4428  !e3_hat(3) = g3(3)/g3_abs
4429 
4430  !--------------------------------------------------
4431 
4432  ! Jacobian
4433  !det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
4434  ! +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
4435  ! +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )
4436 
4437  !--------------------------------------------------
4438 
4439  g1_cross_g2(1) = g1(2)*g2(3)-g1(3)*g2(2)
4440  g1_cross_g2(2) = g1(3)*g2(1)-g1(1)*g2(3)
4441  g1_cross_g2(3) = g1(1)*g2(2)-g1(2)*g2(1)
4442 
4443  !--------------------------------------------------
4444 
4445  do nb = 1, nn
4446 
4447  jsize1 = ndof*(nb-1)+1
4448  jsize2 = ndof*(nb-1)+2
4449  jsize3 = ndof*(nb-1)+3
4450  jsize4 = ndof*(nb-1)+4
4451  jsize5 = ndof*(nb-1)+5
4452  jsize6 = ndof*(nb-1)+6
4453 
4454  n(1, jsize1) = shapefunc(nb)
4455  n(1, jsize2) = 0.0d0
4456  n(1, jsize3) = 0.0d0
4457  n(1, jsize4) = 0.0d0
4458  n(1, jsize5) = 0.0d0
4459  n(1, jsize6) = 0.0d0
4460  n(2, jsize1) = 0.0d0
4461  n(2, jsize2) = shapefunc(nb)
4462  n(2, jsize3) = 0.0d0
4463  n(2, jsize4) = 0.0d0
4464  n(2, jsize5) = 0.0d0
4465  n(2, jsize6) = 0.0d0
4466  n(3, jsize1) = 0.0d0
4467  n(3, jsize2) = 0.0d0
4468  n(3, jsize3) = shapefunc(nb)
4469  n(3, jsize4) = 0.0d0
4470  n(3, jsize5) = 0.0d0
4471  n(3, jsize6) = 0.0d0
4472 
4473  end do
4474 
4475  do isize = 1, ndof*nn
4476 
4477  vect(isize) &
4478  = vect(isize) &
4479  +w_w_lx*( n(1, isize)*g1_cross_g2(1) &
4480  +n(2, isize)*g1_cross_g2(2) &
4481  +n(3, isize)*g1_cross_g2(3) )*val
4482 
4483  end do
4484 
4485  !--------------------------------------------------
4486 
4487  end do
4488 
4489  !--------------------------------------------------------
4490 
4491  end if
4492 
4493  !--------------------------------------------------------------------
4494 
4495  !** Volume load
4496  if( ivol .EQ. 1 ) then
4497 
4498  !--------------------------------------------------------
4499  n_totlyr = gausses(1)%pMaterial%totallyr
4500  do n_layer=1,n_totlyr
4501  do ly = 1, ny
4502 
4503  !--------------------------------------------------
4504 
4505 
4506  call fstr_shell_layer_quadrature_gauss( etype, gausses(1), n_layer, ly, &
4507  zeta_ly, w_ly, ierr_quad )
4508  if( ierr_quad /= 0 ) cycle
4509 
4510  !--------------------------------------------------
4511 
4512  do lx = 1, numofquadpoints(fetype)
4513 
4514  !--------------------------------------------
4515 
4516  call getquadpoint(fetype, lx, naturalcoord)
4517 
4518  xi_lx = naturalcoord(1)
4519  eta_lx = naturalcoord(2)
4520 
4521  w_w_lx = getweight(fetype, lx)
4522 
4523  call getshapefunc(fetype, naturalcoord, shapefunc)
4524 
4525  call getshapederiv(fetype, naturalcoord, shapederiv)
4526 
4527  !--------------------------------------------
4528 
4529  do na = 1, nn
4530 
4531  do i = 1, 3
4532 
4533  u_rot(i, na) &
4534  = shapefunc(na) &
4535  *( zeta_ly*a_over_2_v3(i, na) )
4536 
4537  dudxi_rot(i, na) &
4538  = shapederiv(na, 1) &
4539  *( zeta_ly*a_over_2_v3(i, na) )
4540  dudeta_rot(i, na) &
4541  = shapederiv(na, 2) &
4542  *( zeta_ly*a_over_2_v3(i, na) )
4543  dudzeta_rot(i, na) &
4544  = shapefunc(na) &
4545  *( a_over_2_v3(i, na) )
4546 
4547  end do
4548 
4549  end do
4550 
4551  !--------------------------------------------
4552 
4553  ! Covariant basis vector
4554  do i = 1, 3
4555 
4556  g1(i) = 0.0d0
4557  g2(i) = 0.0d0
4558  g3(i) = 0.0d0
4559 
4560  do na = 1, nn
4561 
4562  g1(i) = g1(i)+shapederiv(na, 1) &
4563  *elem(i, na) &
4564  +dudxi_rot(i, na)
4565  g2(i) = g2(i)+shapederiv(na, 2) &
4566  *elem(i, na) &
4567  +dudeta_rot(i, na)
4568  g3(i) = g3(i)+dudzeta_rot(i, na)
4569 
4570  end do
4571 
4572  end do
4573 
4574  !--------------------------------------------
4575 
4576  ! Jacobian
4577  det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
4578  +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
4579  +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )
4580 
4581  !--------------------------------------------
4582 
4583  ! [ N ] matrix
4584  do nb = 1, nn
4585 
4586  jsize1 = ndof*(nb-1)+1
4587  jsize2 = ndof*(nb-1)+2
4588  jsize3 = ndof*(nb-1)+3
4589  jsize4 = ndof*(nb-1)+4
4590  jsize5 = ndof*(nb-1)+5
4591  jsize6 = ndof*(nb-1)+6
4592 
4593  n(1, jsize1) = shapefunc(nb)
4594  n(2, jsize1) = 0.0d0
4595  n(3, jsize1) = 0.0d0
4596  n(1, jsize2) = 0.0d0
4597  n(2, jsize2) = shapefunc(nb)
4598  n(3, jsize2) = 0.0d0
4599  n(1, jsize3) = 0.0d0
4600  n(2, jsize3) = 0.0d0
4601  n(3, jsize3) = shapefunc(nb)
4602  n(1, jsize4) = 0.0d0
4603  n(2, jsize4) = -u_rot(3, nb)
4604  n(3, jsize4) = u_rot(2, nb)
4605  n(1, jsize5) = u_rot(3, nb)
4606  n(2, jsize5) = 0.0d0
4607  n(3, jsize5) = -u_rot(1, nb)
4608  n(1, jsize6) = -u_rot(2, nb)
4609  n(2, jsize6) = u_rot(1, nb)
4610  n(3, jsize6) = 0.0d0
4611 
4612  enddo
4613 
4614  !--------------------------------------------
4615 
4616  w_w_w_det = w_w_lx*w_ly*det
4617 
4618  !--------------------------------------------
4619 
4620  if( ltype .EQ. 1 ) then
4621 
4622  do isize = 1, ndof*nn
4623 
4624  vect(isize) = vect(isize)+w_w_w_det*n(1, isize)*val
4625 
4626  end do
4627 
4628  else if( ltype .EQ. 2 ) then
4629 
4630  do isize = 1, ndof*nn
4631 
4632  vect(isize) = vect(isize)+w_w_w_det*n(2, isize)*val
4633 
4634  end do
4635 
4636  else if( ltype .EQ. 3 ) then
4637 
4638  do isize = 1, ndof*nn
4639 
4640  vect(isize) = vect(isize)+w_w_w_det*n(3, isize)*val
4641 
4642  end do
4643 
4644  else if( ltype .EQ. 4 ) then
4645 
4646  do isize = 1, ndof*nn
4647 
4648  vect(isize) = vect(isize)+w_w_w_det*rho*ax*n(1, isize)*val
4649  vect(isize) = vect(isize)+w_w_w_det*rho*ay*n(2, isize)*val
4650  vect(isize) = vect(isize)+w_w_w_det*rho*az*n(3, isize)*val
4651 
4652  end do
4653 
4654  else if( ltype .EQ. 5 ) then
4655 
4656  x = 0.0d0
4657  y = 0.0d0
4658  z = 0.0d0
4659 
4660  do nb = 1, nn
4661 
4662  x = x+shapefunc(nb)*elem(1, nb)
4663  y = y+shapefunc(nb)*elem(2, nb)
4664  z = z+shapefunc(nb)*elem(3, nb)
4665 
4666  end do
4667 
4668  hx = ax+( (x-ax)*rx+(y-ay)*ry+(z-az)*rz )/( rx**2+ry**2+rz**2 )*rx
4669  hy = ay+( (x-ax)*rx+(y-ay)*ry+(z-az)*rz )/( rx**2+ry**2+rz**2 )*ry
4670  hz = az+( (x-ax)*rx+(y-ay)*ry+(z-az)*rz )/( rx**2+ry**2+rz**2 )*rz
4671 
4672  phx = x-hx
4673  phy = y-hy
4674  phz = z-hz
4675 
4676  coefx = phx*val*rho*val
4677  coefy = phy*val*rho*val
4678  coefz = phz*val*rho*val
4679 
4680  do isize = 1, ndof*nn
4681 
4682  vect(isize) &
4683  = vect(isize) &
4684  +w_w_w_det*( n(1, isize)*coefx &
4685  +n(2, isize)*coefy &
4686  +n(3, isize)*coefz )
4687 
4688  end do
4689 
4690  end if
4691 
4692  !--------------------------------------------
4693 
4694  end do
4695 
4696  !----------------------------------------------
4697 
4698  end do
4699 
4700  !----------------------------------------------
4701 
4702  end do
4703 
4704  !--------------------------------------------------------
4705 
4706  end if
4707  !--------------------------------------------------------------------
4708 
4709  return
4710 
4711  !####################################################################
4712  end subroutine dl_shell
4713  !####################################################################
4714 
4715 
4716  !####################################################################
4717  subroutine dl_shell_33 &
4718  (ic_type, nn, ndof, xx, yy, zz, rho, thick, &
4719  ltype, params, vect, nsize, gausses)
4720  !####################################################################
4721 
4722  use hecmw
4723  use m_utilities
4724  use mmechgauss
4725 
4726  type(tgaussstatus) :: gausses(:)
4727  !--------------------------------------------------------------------
4728 
4729  integer(kind = kint) :: ic_type
4730  integer(kind = kint) :: nn
4731  integer(kind = kint) :: ndof
4732  real(kind = kreal) :: xx(*), yy(*), zz(*)
4733  real(kind = kreal) :: rho
4734  real(kind = kreal) :: thick
4735  real(kind = kreal) :: params(*)
4736  real(kind = kreal) :: vect(*)
4737  integer(kind = kint) :: nsize
4738  integer :: ltype, i
4739  real(kind = kreal) :: tmp(24)
4740  !--------------------------------------------------------------------
4741 
4742  if(ic_type == 761)then
4743  !ic_type = 731
4744  !nn = 3
4745  !ndof = 6
4746  call dl_shell(731, 3, 6, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
4747  !ic_type = 761
4748  !nn = 6
4749  !ndof = 3
4750 
4751  tmp = 0.0
4752  do i=1,18
4753  tmp(i) = vect(i)
4754  enddo
4755 
4756  vect( 1) = tmp(1)
4757  vect( 2) = tmp(2)
4758  vect( 3) = tmp(3)
4759  vect( 4) = tmp(7)
4760  vect( 5) = tmp(8)
4761  vect( 6) = tmp(9)
4762  vect( 7) = tmp(13)
4763  vect( 8) = tmp(14)
4764  vect( 9) = tmp(15)
4765  vect(10) = tmp(4)
4766  vect(11) = tmp(5)
4767  vect(12) = tmp(6)
4768  vect(13) = tmp(10)
4769  vect(14) = tmp(11)
4770  vect(15) = tmp(12)
4771  vect(16) = tmp(16)
4772  vect(17) = tmp(17)
4773  vect(18) = tmp(18)
4774 
4775  elseif(ic_type == 781)then
4776  !ic_type = 741
4777  !nn = 4
4778  !ndof = 6
4779  call dl_shell(741, 4, 6, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
4780  !ic_type = 781
4781  !nn = 8
4782  !ndof = 3
4783 
4784  tmp = 0.0
4785  do i=1,24
4786  tmp(i) = vect(i)
4787  enddo
4788 
4789  vect( 1) = tmp(1)
4790  vect( 2) = tmp(2)
4791  vect( 3) = tmp(3)
4792  vect( 4) = tmp(7)
4793  vect( 5) = tmp(8)
4794  vect( 6) = tmp(9)
4795  vect( 7) = tmp(13)
4796  vect( 8) = tmp(14)
4797  vect( 9) = tmp(15)
4798  vect(10) = tmp(19)
4799  vect(11) = tmp(20)
4800  vect(12) = tmp(21)
4801  vect(13) = tmp(4)
4802  vect(14) = tmp(5)
4803  vect(15) = tmp(6)
4804  vect(16) = tmp(10)
4805  vect(17) = tmp(11)
4806  vect(18) = tmp(12)
4807  vect(19) = tmp(16)
4808  vect(20) = tmp(17)
4809  vect(21) = tmp(18)
4810  vect(22) = tmp(22)
4811  vect(23) = tmp(23)
4812  vect(24) = tmp(24)
4813 
4814  endif
4815 
4816  end subroutine dl_shell_33
4817 
4818  !####################################################################
4819  ! Update shell stress and equivalent nodal force.
4820  subroutine updatest_shell_mitc &
4821  (etype, nn, ndof, ecoord, u, du, gausses, qf, thick, mixflag, nddisp, element, nddirector, ndrefdirector, &
4822  ndcurdirector, nddrill)
4823  !####################################################################
4824 
4825  use mmechgauss
4826  use m_matmatrix
4827  use mmaterial, only: totallag, updatelag, iselastic
4828 
4829  !--------------------------------------------------------------------
4830 
4831  integer(kind = kint), intent(in) :: etype
4832  integer(kind = kint), intent(in) :: nn, mixflag
4833  integer(kind = kint), intent(in) :: ndof
4834  real(kind = kreal), intent(in) :: ecoord(3, nn)
4835  real(kind = kreal), intent(in) :: u(:, :)
4836  real(kind = kreal), intent(in) :: du(:, :)
4837  type(tgaussstatus), intent(in) :: gausses(:)
4838  real(kind = kreal), intent(out) :: qf(:)
4839  real(kind = kreal), intent(in) :: thick
4840 
4841  real(kind = kreal), intent(in), optional :: nddisp(ndof, nn)
4842  type(telement), intent(inout), optional :: element
4843  real(kind = kreal), intent(in), optional :: nddirector(3, nn)
4844  real(kind = kreal), intent(in), optional :: ndrefdirector(3, nn)
4845  real(kind = kreal), intent(in), optional :: ndcurdirector(3, nn)
4846  real(kind = kreal), intent(in), optional :: nddrill(nn)
4847  !--------------------------------------------------------------------
4848 
4849  real(kind = kreal) :: stiff(nn*ndof, nn*ndof), totaldisp(nn*ndof), edisp(6, nn), qf_direct(nn*ndof)
4850  real(kind = kreal) :: incdisp(6, nn), basedisp(6, nn)
4851  integer(kind = kint) :: i
4852  integer(kind = kint) :: flag
4853  logical :: use_stress_force
4854 
4855  flag = gausses(1)%pMaterial%nlgeom_flag
4856  use_stress_force = present( element ) .and. &
4857  shellsupportsfiniterotationkinematics( etype, nn ) .and. &
4858  ( flag == totallag .or. flag == updatelag ) .and. &
4859  iselastic( gausses(1)%pMaterial%mtype )
4860 
4861  totaldisp = 0.d0
4862  edisp(:, :) = 0.0d0
4863  incdisp(:, :) = 0.0d0
4864  basedisp(:, :) = 0.0d0
4865  do i=1,nn
4866  totaldisp(ndof*(i-1)+1:ndof*i) = u(1:ndof,i) + du(1:ndof,i)
4867  basedisp(1:6, i) = u(1:6, i)
4868  incdisp(1:6, i) = du(1:6, i)
4869  if( use_stress_force ) then
4870  edisp(1:3, i) = u(1:3, i) + du(1:3, i)
4871  call shellcomposerotationvector( u(4:6, i), du(4:6, i), edisp(4:6, i) )
4872  else
4873  edisp(1:6, i) = u(1:6, i) + du(1:6, i)
4874  endif
4875  end do
4876 
4877  ! Equivalent nodal force from layer stress.
4878  qf_direct(:) = 0.0d0
4879 
4880  if( use_stress_force ) then
4881  if( flag == updatelag ) then
4882  if( present( ndcurdirector ) ) then
4883  call updateshelllayergauss_shell_mitc( &
4884  etype, nn, ndof, ecoord, element, incdisp, thick, nddirector=nddirector, &
4885  ndrefdirector=ndcurdirector, ndbase_disp=basedisp )
4886  else
4887  call updateshelllayergauss_shell_mitc( &
4888  etype, nn, ndof, ecoord, element, incdisp, thick, nddirector=nddirector, &
4889  ndrefdirector=ndrefdirector, ndbase_disp=basedisp )
4890  endif
4891  else
4892  call updateshelllayergauss_shell_mitc( &
4893  etype, nn, ndof, ecoord, element, edisp, thick, nddirector=nddirector, ndrefdirector=ndrefdirector )
4894  endif
4895  if( present( nddisp ) ) then
4896  call stf_shell_mitc(etype, nn, ndof, ecoord, gausses, stiff, thick, &
4897  mixflag, nddisp=nddisp, element=element, qf_stress=qf_direct, &
4898  include_geo_stiff=.false., nddirector=nddirector, ndrefdirector=ndrefdirector, &
4899  ndcurdirector=ndcurdirector, nddrill=nddrill)
4900  else
4901  call stf_shell_mitc(etype, nn, ndof, ecoord, gausses, stiff, thick, &
4902  mixflag, nddisp=edisp, element=element, qf_stress=qf_direct, &
4903  include_geo_stiff=.false., nddirector=nddirector, ndrefdirector=ndrefdirector, &
4904  ndcurdirector=ndcurdirector, nddrill=nddrill)
4905  endif
4906  else
4907  if( present( nddisp ) ) then
4908  call stf_shell_mitc(etype, nn, ndof, ecoord, gausses, stiff, thick, &
4909  mixflag, nddisp=nddisp, include_geo_stiff=.false., nddirector=nddirector, ndrefdirector=ndrefdirector, &
4910  ndcurdirector=ndcurdirector, nddrill=nddrill)
4911  else
4912  call stf_shell_mitc(etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag)
4913  endif
4914  endif
4915 
4916  if( use_stress_force ) then
4917  qf = qf_direct
4918  else
4919  qf = matmul(stiff,totaldisp)
4920  endif
4921 
4922  end subroutine updatest_shell_mitc
4923 
4924  !####################################################################
4925  ! this subroutine can be used only for linear analysis
4926  subroutine updatest_shell_mitc33 &
4927  (etype, nn, ndof, ecoord, u, du, gausses, qf, thick, mixflag, nddisp)
4928  !####################################################################
4929 
4930  use mmechgauss
4931  use m_matmatrix
4932 
4933  !--------------------------------------------------------------------
4934 
4935  integer(kind = kint), intent(in) :: etype
4936  integer(kind = kint), intent(in) :: nn, mixflag
4937  integer(kind = kint), intent(in) :: ndof
4938  real(kind = kreal), intent(in) :: ecoord(3, nn)
4939  real(kind = kreal), intent(in) :: u(3, nn*2)
4940  real(kind = kreal), intent(in) :: du(3, nn*2)
4941  type(tgaussstatus), intent(in) :: gausses(:)
4942  real(kind = kreal), intent(out) :: qf(:)
4943  real(kind = kreal), intent(in) :: thick
4944 
4945  real(kind = kreal), intent(in), optional :: nddisp(3, nn)
4946  !--------------------------------------------------------------------
4947 
4948  real(kind = kreal) :: stiff(nn*ndof, nn*ndof), totaldisp(nn*ndof)
4949  integer(kind = kint) :: i
4950 
4951  call stf_shell_mitc(etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag, nddisp)
4952 
4953  totaldisp = 0.d0
4954  do i=1,nn
4955  totaldisp(ndof*(i-1)+1:ndof*(i-1)+3) = u(1:3,2*i-1) + du(1:3,2*i-1)
4956  totaldisp(ndof*(i-1)+4:ndof*(i-1)+6) = u(1:3,2*i) + du(1:3,2*i)
4957  end do
4958 
4959  qf = matmul(stiff,totaldisp)
4960 
4961  end subroutine updatest_shell_mitc33
4962 
4963  !####################################################################
4964  subroutine mass_shell(etype, nn, elem, rho, thick, gausses, mass, lumped)
4965  !####################################################################
4966  use hecmw
4967  use m_utilities
4968  use mmechgauss
4969  use quadrature
4970  integer(kind = kint), intent(in) :: etype
4971  integer(kind = kint), intent(in) :: nn
4972  real(kind = kreal), intent(in) :: elem(3,nn)
4973  real(kind = kreal), intent(in) :: rho
4974  real(kind = kreal), intent(in) :: thick
4975  type(tgaussstatus), intent(in) :: gausses(:)
4976  real(kind=kreal), intent(out) :: mass(:,:)
4977  real(kind=kreal), intent(out) :: lumped(:)
4978 
4979  !--------------------------------------------------------------------
4980 
4981  integer :: lx, ly, nsize, ndof
4982  integer :: fetype
4983  integer :: ny
4984  integer :: i
4985  integer(kind=kint) :: ierr_quad
4986  integer :: na, nb
4987  integer :: jsize1, jsize2, jsize3, jsize4, jsize5, jsize6
4988  integer :: n_totlyr, n_layer
4989 
4990  real(kind = kreal) :: xi_lx, eta_lx, zeta_ly
4991  real(kind = kreal) :: w_w_lx, w_ly
4992  real(kind = kreal) :: naturalcoord(2)
4993  real(kind = kreal) :: nncoord(nn, 2)
4994  real(kind = kreal) :: shapefunc(nn)
4995  real(kind = kreal) :: shapederiv(nn, 2)
4996  real(kind = kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
4997  real(kind = kreal) :: v1_abs, v2_abs, v3_abs
4998  real(kind = kreal) :: a_over_2_v3(3, nn)
4999  real(kind = kreal) :: u_rot(3, nn)
5000  real(kind = kreal) :: dudxi_rot(3, nn), dudeta_rot(3, nn), dudzeta_rot(3, nn)
5001  real(kind = kreal) :: g1(3), g2(3), g3(3)
5002  real(kind = kreal) :: e_0(3)
5003  real(kind = kreal) :: det
5004  real(kind = kreal) :: det_cg3(3)
5005  real(kind = kreal) :: det_cg3_abs
5006  real(kind = kreal) :: w_w_w_det
5007  real(kind = kreal) :: n(3, 6*nn)
5008  real(kind = kreal) :: totalmass, totdiag
5009 
5010  ny = 0; ndof=6
5011  nsize = ndof*nn
5012 
5013  !--------------------------------------------------------------------
5014 
5015  ! MITC4
5016  if( etype == fe_mitc4_shell ) then
5017  fetype = fe_mitc4_shell
5018  ny = 2
5019 
5020  ! MITC9
5021  else if( etype == fe_mitc9_shell ) then
5022  fetype = fe_mitc9_shell
5023  ny = 3
5024 
5025  ! MITC3
5026  else if( etype == fe_mitc3_shell ) then
5027  fetype = fe_mitc3_shell
5028  ny = 2
5029 
5030  end if
5031 
5032  !-------------------------------------------------------------------
5033 
5034  ! xi-coordinate at a node in a local element
5035  ! eta-coordinate at a node in a local element
5036  call getnodalnaturalcoord(fetype, nncoord)
5037 
5038  ! xi-coordinate at the center point in a local element
5039  ! eta-coordinate at the center point in a local element
5040  naturalcoord(1) = 0.0d0
5041  naturalcoord(2) = 0.0d0
5042 
5043  call getshapederiv(fetype, naturalcoord, shapederiv)
5044 
5045  !--------------------------------------------------------------
5046 
5047  ! Covariant basis vector
5048  g1(:) = matmul( elem, shapederiv(:,1) )
5049 
5050  e_0(:) = g1(:)
5051 
5052  !--------------------------------------------------------------
5053 
5054  do nb = 1, nn
5055 
5056  !--------------------------------------------------------
5057 
5058  naturalcoord(1) = nncoord(nb, 1)
5059  naturalcoord(2) = nncoord(nb, 2)
5060  call getshapederiv(fetype, naturalcoord, shapederiv)
5061 
5062  !--------------------------------------------------------
5063 
5064  ! Covariant basis vector
5065  g1(:) = matmul( elem, shapederiv(:,1) )
5066  g2(:) = matmul( elem, shapederiv(:,2) )
5067 
5068  !--------------------------------------------------------
5069 
5070  det_cg3(1) = g1(2)*g2(3)-g1(3)*g2(2)
5071  det_cg3(2) = g1(3)*g2(1)-g1(1)*g2(3)
5072  det_cg3(3) = g1(1)*g2(2)-g1(2)*g2(1)
5073 
5074  det_cg3_abs = dsqrt( det_cg3(1)*det_cg3(1) &
5075  +det_cg3(2)*det_cg3(2) &
5076  +det_cg3(3)*det_cg3(3) )
5077 
5078  v3(:, nb) = det_cg3(:)/det_cg3_abs
5079 
5080  !--------------------------------------------------------
5081 
5082  v2(1, nb) = v3(2, nb)*e_0(3)-v3(3, nb)*e_0(2)
5083  v2(2, nb) = v3(3, nb)*e_0(1)-v3(1, nb)*e_0(3)
5084  v2(3, nb) = v3(1, nb)*e_0(2)-v3(2, nb)*e_0(1)
5085 
5086  v2_abs = dsqrt( v2(1, nb)*v2(1, nb) &
5087  +v2(2, nb)*v2(2, nb) &
5088  +v2(3, nb)*v2(3, nb) )
5089 
5090  if( v2_abs > 1.0d-15 ) then
5091 
5092  v2(1, nb) = v2(1, nb)/v2_abs
5093  v2(2, nb) = v2(2, nb)/v2_abs
5094  v2(3, nb) = v2(3, nb)/v2_abs
5095 
5096  v1(1, nb) = v2(2, nb)*v3(3, nb) &
5097  -v2(3, nb)*v3(2, nb)
5098  v1(2, nb) = v2(3, nb)*v3(1, nb) &
5099  -v2(1, nb)*v3(3, nb)
5100  v1(3, nb) = v2(1, nb)*v3(2, nb) &
5101  -v2(2, nb)*v3(1, nb)
5102 
5103  v1_abs = dsqrt( v1(1, nb)*v1(1, nb) &
5104  +v1(2, nb)*v1(2, nb) &
5105  +v1(3, nb)*v1(3, nb) )
5106 
5107  v1(1, nb) = v1(1, nb)/v1_abs
5108  v1(2, nb) = v1(2, nb)/v1_abs
5109  v1(3, nb) = v1(3, nb)/v1_abs
5110 
5111  else
5112 
5113  v1(1, nb) = 0.0d0
5114  v1(2, nb) = 0.0d0
5115  v1(3, nb) = -1.0d0
5116 
5117  v2(1, nb) = 0.0d0
5118  v2(2, nb) = 1.0d0
5119  v2(3, nb) = 0.0d0
5120 
5121  end if
5122 
5123  !--------------------------------------------------------
5124 
5125  v3(1, nb) = v1(2, nb)*v2(3, nb) &
5126  -v1(3, nb)*v2(2, nb)
5127  v3(2, nb) = v1(3, nb)*v2(1, nb) &
5128  -v1(1, nb)*v2(3, nb)
5129  v3(3, nb) = v1(1, nb)*v2(2, nb) &
5130  -v1(2, nb)*v2(1, nb)
5131 
5132  v3_abs = dsqrt( v3(1, nb)*v3(1, nb) &
5133  +v3(2, nb)*v3(2, nb) &
5134  +v3(3, nb)*v3(3, nb) )
5135 
5136  v3(1, nb) = v3(1, nb)/v3_abs
5137  v3(2, nb) = v3(2, nb)/v3_abs
5138  v3(3, nb) = v3(3, nb)/v3_abs
5139 
5140  !--------------------------------------------------------
5141 
5142  a_over_2_v3(1, nb) = 0.5d0*thick*v3(1, nb)
5143  a_over_2_v3(2, nb) = 0.5d0*thick*v3(2, nb)
5144  a_over_2_v3(3, nb) = 0.5d0*thick*v3(3, nb)
5145 
5146  !--------------------------------------------------------
5147 
5148  end do
5149 
5150  !--------------------------------------------------------------------
5151 
5152  mass(:,:) = 0.0d0
5153  totalmass = 0.d0
5154  n_totlyr = gausses(1)%pMaterial%totallyr
5155  do n_layer=1,n_totlyr
5156  do ly = 1, ny
5157 
5158  !--------------------------------------------------
5159 
5160  call fstr_shell_layer_quadrature_gauss( etype, gausses(1), n_layer, ly, &
5161  zeta_ly, w_ly, ierr_quad )
5162  if( ierr_quad /= 0 ) cycle
5163 
5164  !--------------------------------------------------
5165 
5166  do lx = 1, numofquadpoints(fetype)
5167 
5168  !--------------------------------------------
5169 
5170  call getquadpoint(fetype, lx, naturalcoord)
5171 
5172  xi_lx = naturalcoord(1)
5173  eta_lx = naturalcoord(2)
5174 
5175  w_w_lx = getweight(fetype, lx)
5176 
5177  call getshapefunc(fetype, naturalcoord, shapefunc)
5178  call getshapederiv(fetype, naturalcoord, shapederiv)
5179 
5180  !--------------------------------------------
5181 
5182  do na = 1, nn
5183 
5184  do i = 1, 3
5185 
5186  u_rot(i, na) = shapefunc(na)*( zeta_ly*a_over_2_v3(i, na) )
5187 
5188  dudxi_rot(i, na) = shapederiv(na, 1) &
5189  *( zeta_ly*a_over_2_v3(i, na) )
5190  dudeta_rot(i, na) = shapederiv(na, 2) &
5191  *( zeta_ly*a_over_2_v3(i, na) )
5192  dudzeta_rot(i, na) = shapefunc(na) &
5193  *( a_over_2_v3(i, na) )
5194 
5195  end do
5196 
5197  end do
5198 
5199  !--------------------------------------------
5200 
5201  ! Covariant basis vector
5202  do i = 1, 3
5203  g1(i) = 0.0d0
5204  g2(i) = 0.0d0
5205  g3(i) = 0.0d0
5206  do na = 1, nn
5207  g1(i) = g1(i)+shapederiv(na, 1) *elem(i, na) &
5208  +dudxi_rot(i, na)
5209  g2(i) = g2(i)+shapederiv(na, 2) *elem(i, na) &
5210  +dudeta_rot(i, na)
5211  g3(i) = g3(i)+dudzeta_rot(i, na)
5212  end do
5213  end do
5214 
5215  !--------------------------------------------
5216 
5217  ! Jacobian
5218  det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
5219  +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
5220  +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )
5221 
5222  !--------------------------------------------
5223 
5224  ! [ N ] matrix
5225  do nb = 1, nn
5226 
5227  jsize1 = ndof*(nb-1)+1
5228  jsize2 = ndof*(nb-1)+2
5229  jsize3 = ndof*(nb-1)+3
5230  jsize4 = ndof*(nb-1)+4
5231  jsize5 = ndof*(nb-1)+5
5232  jsize6 = ndof*(nb-1)+6
5233 
5234  n(1, jsize1) = shapefunc(nb)
5235  n(2, jsize1) = 0.0d0
5236  n(3, jsize1) = 0.0d0
5237  n(1, jsize2) = 0.0d0
5238  n(2, jsize2) = shapefunc(nb)
5239  n(3, jsize2) = 0.0d0
5240  n(1, jsize3) = 0.0d0
5241  n(2, jsize3) = 0.0d0
5242  n(3, jsize3) = shapefunc(nb)
5243  n(1, jsize4) = 0.0d0
5244  n(2, jsize4) = -u_rot(3, nb)
5245  n(3, jsize4) = u_rot(2, nb)
5246  n(1, jsize5) = u_rot(3, nb)
5247  n(2, jsize5) = 0.0d0
5248  n(3, jsize5) = -u_rot(1, nb)
5249  n(1, jsize6) = -u_rot(2, nb)
5250  n(2, jsize6) = u_rot(1, nb)
5251  n(3, jsize6) = 0.0d0
5252 
5253  enddo
5254 
5255  !--------------------------------------------
5256 
5257  w_w_w_det = w_w_lx*w_ly*det*gausses(1)%pMaterial%shell_var(n_layer)%weight
5258  mass(1:nsize,1:nsize) = mass(1:nsize,1:nsize)+ matmul( transpose(n), n )*w_w_w_det*rho
5259  totalmass = totalmass + w_w_w_det*rho
5260  !--------------------------------------------
5261 
5262  end do
5263 
5264  !----------------------------------------------
5265 
5266  end do
5267 
5268  !----------------------------------------------
5269 
5270  end do
5271  totalmass = totalmass*3.d0
5272 
5273  totdiag=0.d0
5274  do nb = 1, nn
5275  DO i = 1, 3
5276  lx = (nb-1)*ndof+i
5277  totdiag = totdiag + mass(lx,lx)
5278  END DO
5279  ENDDO
5280  do nb = 1, nn
5281  DO i = 1, 6
5282  lx = (nb-1)*ndof+i
5283  lumped(lx) = mass(lx,lx)/totdiag* totalmass
5284  END DO
5285  ENDDO
5286 
5287  !####################################################################
5288  end subroutine mass_shell
5289  !####################################################################
5290 
5291 end module m_static_lib_shell
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
elementinfo::getshapefunc
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coordinate system.
Definition: element.f90:647
elementinfo::getquadpoint
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:489
elementinfo::getshapederiv
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate derivatives of shape function in natural coordinate system.
Definition: element.f90:578
elementinfo::fe_mitc4_shell
integer, parameter fe_mitc4_shell
Definition: element.f90:94
elementinfo::fe_mitc9_shell
integer, parameter fe_mitc9_shell
Definition: element.f90:96
elementinfo::getweight
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:535
elementinfo::numofquadpoints
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:450
elementinfo::fe_mitc3_shell
integer, parameter fe_mitc3_shell
Definition: element.f90:93
elementinfo::getnodalnaturalcoord
subroutine getnodalnaturalcoord(fetype, nncoord)
Definition: element.f90:699
hecmw
Definition: hecmw.f90:6
m_fstr_finiterotationkinematics
Shared predicates for finite-rotation nodal kinematics.
Definition: fstr_FiniteRotationKinematics.f90:7
m_fstr_finiterotationkinematics::fstr_is_finite_rotation_shell_element
logical function, public fstr_is_finite_rotation_shell_element(etype, nn)
Definition: fstr_FiniteRotationKinematics.f90:24
m_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.F90:15
m_fstr::kopss_solution
integer(kind=kint), parameter kopss_solution
Definition: m_fstr.F90:133
m_fstr::opsstype
integer(kind=kint) opsstype
Definition: m_fstr.F90:135
m_matmatrix
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
m_matmatrix::matlmatrix_shell
subroutine matlmatrix_shell(gauss, sectType, D, e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, alpha, n_layer)
Definition: calMatMatrix.f90:220
m_static_lib_shell
This module ...
Definition: static_LIB_shell.f90:6
m_static_lib_shell::stf_shell_mitc
subroutine stf_shell_mitc(etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag, nddisp, element, qf_stress, include_geo_stiff, nddirector, ndrefdirector, ndcurdirector, nddrill)
Definition: static_LIB_shell.f90:285
m_static_lib_shell::outer_product3
pure real(kind=kreal) function, dimension(3, 3) outer_product3(a, b)
Definition: static_LIB_shell.f90:64
m_static_lib_shell::shellstressvectortotensor
subroutine shellstressvectortotensor(stress, tensor)
Definition: static_LIB_shell.f90:118
m_static_lib_shell::shelldirectorincrement
subroutine shelldirectorincrement(theta, director_ref, director_inc)
Definition: static_LIB_shell.f90:2273
m_static_lib_shell::shelldirectorincrementalsecondderiv
pure subroutine shelldirectorincrementalsecondderiv(director_current, director_second)
Definition: static_LIB_shell.f90:2600
m_static_lib_shell::elementstress_shell_mitc
subroutine elementstress_shell_mitc(etype, nn, ndof, ecoord, gausses, edisp, strain, stress, thick, zeta, n_layer, n_totlyr, surface_gauss_points, local_strain, local_stress, local_stress_override, nddirector, ndrefdirector, ndbase_disp)
Definition: static_LIB_shell.f90:2646
m_static_lib_shell::shelltensortostressvector
subroutine shelltensortostressvector(tensor, stress)
Definition: static_LIB_shell.f90:137
m_static_lib_shell::shell_basis_from_covariant
pure subroutine shell_basis_from_covariant(g1, g2, g3, e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, det)
Definition: static_LIB_shell.f90:79
m_static_lib_shell::updatestressshellul_elastic
subroutine updatestressshellul_elastic(gauss, dstrain, dstress, trace_coeff)
Definition: static_LIB_shell.f90:3928
m_static_lib_shell::shelladdulobjectivetracetangent
subroutine shelladdulobjectivetracetangent(stress_old, ncol, B, DB, trace_coeff)
Definition: static_LIB_shell.f90:245
m_static_lib_shell::shellrelativerotationvector
pure subroutine shellrelativerotationvector(theta_old, theta_target, theta_inc)
Definition: static_LIB_shell.f90:2425
m_static_lib_shell::dl_shell_33
subroutine dl_shell_33(ic_type, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
Definition: static_LIB_shell.f90:4720
m_static_lib_shell::shellobjectivestressincrement
subroutine shellobjectivestressincrement(stress_old, dstrain, dstress_obj, trace_coeff)
Definition: static_LIB_shell.f90:153
m_static_lib_shell::shellorthonormalizetriad
pure subroutine shellorthonormalizetriad(triad_in, triad_out)
Definition: static_LIB_shell.f90:2447
m_static_lib_shell::shellcomposenodaldisplacement
subroutine shellcomposenodaldisplacement(ndof, nn, disp_old, disp_inc, disp_new)
Definition: static_LIB_shell.f90:2542
m_static_lib_shell::shellcomposerotationvector
pure subroutine shellcomposerotationvector(theta_old, theta_inc, theta_new)
Definition: static_LIB_shell.f90:2403
m_static_lib_shell::shellrotationvectortomatrix
pure subroutine shellrotationvectortomatrix(theta, rotmat)
Definition: static_LIB_shell.f90:2311
m_static_lib_shell::shellplanestresstracecoeff
real(kind=kreal) function shellplanestresstracecoeff(gauss, n_layer)
Definition: static_LIB_shell.f90:185
m_static_lib_shell::updatest_shell_mitc33
subroutine updatest_shell_mitc33(etype, nn, ndof, ecoord, u, du, gausses, qf, thick, mixflag, nddisp)
Definition: static_LIB_shell.f90:4928
m_static_lib_shell::shellrotationmatrixtovector
pure subroutine shellrotationmatrixtovector(rotmat, theta)
Definition: static_LIB_shell.f90:2356
m_static_lib_shell::updateshelllayergauss_shell_mitc
subroutine updateshelllayergauss_shell_mitc(etype, nn, ndof, ecoord, element, edisp, thick, nddirector, ndrefdirector, ndbase_disp)
Definition: static_LIB_shell.f90:3966
m_static_lib_shell::updatest_shell_mitc
subroutine updatest_shell_mitc(etype, nn, ndof, ecoord, u, du, gausses, qf, thick, mixflag, nddisp, element, nddirector, ndrefdirector, ndcurdirector, nddrill)
Definition: static_LIB_shell.f90:4823
m_static_lib_shell::shelladdulobjectivetangent
subroutine shelladdulobjectivetangent(stress_old, ncol, B, DB, trace_coeff)
Definition: static_LIB_shell.f90:268
m_static_lib_shell::shelladdstressvectortodb
subroutine shelladdstressvectortodb(dstress, j, DB)
Definition: static_LIB_shell.f90:229
m_static_lib_shell::shellsupportsfiniterotationkinematics
logical function shellsupportsfiniterotationkinematics(etype, nn)
Definition: static_LIB_shell.f90:57
m_static_lib_shell::shellobjectivetracestressincrement
subroutine shellobjectivetracestressincrement(stress_old, dstrain, dstress_trace, trace_coeff)
Definition: static_LIB_shell.f90:165
m_static_lib_shell::shellupdatetriadwithincrement
pure subroutine shellupdatetriadwithincrement(triad_old, drill_old, theta_inc, triad_new, drill_new)
Definition: static_LIB_shell.f90:2494
m_static_lib_shell::shelldirectorincrementalderiv
pure subroutine shelldirectorincrementalderiv(director_current, director_deriv)
Definition: static_LIB_shell.f90:2574
m_static_lib_shell::mass_shell
subroutine mass_shell(etype, nn, elem, rho, thick, gausses, mass, lumped)
Definition: static_LIB_shell.f90:4965
m_static_lib_shell::dl_shell
subroutine dl_shell(etype, nn, ndof, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
Definition: static_LIB_shell.f90:4050
m_utilities
This module provides aux functions.
Definition: utilities.f90:6
m_utilities::get_principal
subroutine get_principal(tensor, eigval, princmatrix)
Definition: utilities.f90:450
mmaterial
This module summarizes all information of material properties.
Definition: material.f90:6
mmaterial::getelastictype
integer function getelastictype(mtype)
Get elastic type.
Definition: material.f90:309
mmaterial::totallag
integer(kind=kint), parameter totallag
Definition: material.f90:14
mmaterial::m_poisson
integer(kind=kint), parameter m_poisson
Definition: material.f90:93
mmaterial::infinitesimal
integer(kind=kint), parameter infinitesimal
Definition: material.f90:13
mmaterial::mc_isoelastic
character(len=dict_key_length) mc_isoelastic
Definition: material.f90:140
mmaterial::iselastic
logical function iselastic(mtype)
If it is an elastic material?
Definition: material.f90:354
mmaterial::updatelag
integer(kind=kint), parameter updatelag
Definition: material.f90:15
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
mmechgauss::fstr_shell_layer_zeta
subroutine fstr_shell_layer_zeta(gauss, ilayer, zeta, zeta_layer, ierr)
Map layer-local zeta to the whole shell thickness coordinate.
Definition: mechgauss.f90:227
mmechgauss::fstr_shell_layer_quadrature_gauss
subroutine fstr_shell_layer_quadrature_gauss(etype, gauss, ilayer, ithick, zeta_layer, weight, ierr)
Layer-local shell thickness coordinate and quadrature weight from material status.
Definition: mechgauss.f90:210
mmechgauss::fstr_shell_layer_gauss_index
integer(kind=kint) function fstr_shell_layer_gauss_index(element, ig, ilayer, ithick)
Convert surface Gauss/layer/thickness indices to shell_layer_gausses index.
Definition: mechgauss.f90:139
mmechgauss::fstr_shell_thickness_quadrature
subroutine fstr_shell_thickness_quadrature(etype, ithick, zeta, weight, ierr)
Through-thickness quadrature point and weight used by shell elements.
Definition: mechgauss.f90:160
quadrature
This module contains Gauss point information.
Definition: quadrature.f90:28
quadrature::gauss1d3
real(kind=kreal), dimension(1, 3) gauss1d3
Definition: quadrature.f90:32
mmechgauss::telement
All data should be recorded in every elements.
Definition: mechgauss.f90:35
mmechgauss::tgaussstatus
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:16