FrontISTR  5.9.0
Large-scale structural analysis program with finit element method
static_LIB_beam.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_beam
7 
8  use hecmw
9  use mmechgauss
10  use m_utilities
11  use elementinfo
12  use m_fstr, only: kel611timoshenko
13 
14  implicit none
15 
16 contains
17 
18 
19  subroutine framtr(refx, xl, le, t)
20  real(kind=kreal), intent(in) :: refx(3)
21  real(kind=kreal), intent(in) :: xl(3,2)
22  real(kind=kreal), intent(out) :: le
23  real(kind=kreal), intent(out) :: t(3,3)
24 
25  real(kind=kreal) :: dl
26  real(kind=kreal), parameter :: tol = 1.d-08
27 
28  t(1,1) = xl(1,2) - xl(1,1)
29  t(1,2) = xl(2,2) - xl(2,1)
30  t(1,3) = xl(3,2) - xl(3,1)
31  le = sqrt(t(1,1)*t(1,1)+t(1,2)*t(1,2)+t(1,3)*t(1,3))
32  dl = 1.0d0/le
33  t(1,1) = t(1,1)*dl
34  t(1,2) = t(1,2)*dl
35  t(1,3) = t(1,3)*dl
36 
37  t(3,1) = refx(1)
38  t(3,2) = refx(2)
39  t(3,3) = refx(3)
40 
41 
42  t(2,1) = (t(3,2)*t(1,3) - t(3,3)*t(1,2))
43  t(2,2) = (t(3,3)*t(1,1) - t(3,1)*t(1,3))
44  t(2,3) = (t(3,1)*t(1,2) - t(3,2)*t(1,1))
45  dl = sqrt(t(2,1)*t(2,1)+t(2,2)*t(2,2)+t(2,3)*t(2,3))
46  if(dl<tol*le) then
47  stop "Bad reference for beam element!"
48  else
49  t(2,1) = t(2,1)/dl
50  t(2,2) = t(2,2)/dl
51  t(2,3) = t(2,3)/dl
52  t(3,1) = t(1,2)*t(2,3) - t(1,3)*t(2,2)
53  t(3,2) = t(1,3)*t(2,1) - t(1,1)*t(2,3)
54  t(3,3) = t(1,1)*t(2,2) - t(1,2)*t(2,1)
55  endif
56 
57  end subroutine framtr
58 
60  subroutine stf_beam_local(le, section, E, P, formulation, stiff)
61  real(kind=kreal), intent(in) :: le
62  real(kind=kreal), intent(in) :: section(:)
63  real(kind=kreal), intent(in) :: e, p
64  integer(kind=kint), intent(in), optional :: formulation
65  real(kind=kreal), intent(out) :: stiff(12,12)
66 
67  real(kind=kreal) :: g, l2, l3, a, iy, iz, jx, ea
68  real(kind=kreal) :: phi_y, phi_z
69  real(kind=kreal) :: kyy, kyz, kzz, kryy, kryz
70  real(kind=kreal) :: kww, kwy, krr, krz
71 
72  l2 = le*le
73  l3 = l2*le
74  g = e/(2.d0*(1.d0 + p))
75 
76  a = section(4); iy = section(5); iz = section(6); jx = section(7)
77 
78  phi_y = 0.d0
79  phi_z = 0.d0
80  if( present(formulation) ) then
81  if( formulation == kel611timoshenko ) then
82  ! The BEAM section has no independent shear areas; use A in both planes.
83  phi_y = 12.d0*e*iz/(g*a*l2)
84  phi_z = 12.d0*e*iy/(g*a*l2)
85  endif
86  endif
87 
88  ea = e*a/le
89  kyy = 12.d0*e*iz/(l3*(1.d0 + phi_y))
90  kyz = 6.d0*e*iz/(l2*(1.d0 + phi_y))
91  kzz = (4.d0 + phi_y)*e*iz/(le*(1.d0 + phi_y))
92  kryz = (2.d0 - phi_y)*e*iz/(le*(1.d0 + phi_y))
93  kww = 12.d0*e*iy/(l3*(1.d0 + phi_z))
94  kwy = 6.d0*e*iy/(l2*(1.d0 + phi_z))
95  krr = (4.d0 + phi_z)*e*iy/(le*(1.d0 + phi_z))
96  krz = (2.d0 - phi_z)*e*iy/(le*(1.d0 + phi_z))
97 
98  stiff = 0.d0
99  stiff(1,1) = ea
100  stiff(7,1) = -ea
101 
102  stiff(2,2) = kyy
103  stiff(6,2) = kyz
104  stiff(8,2) = -kyy
105  stiff(12,2) = kyz
106 
107  stiff(3,3) = kww
108  stiff(5,3) = -kwy
109  stiff(9,3) = -kww
110  stiff(11,3) = -kwy
111 
112  stiff(4,4) = g*jx/le
113  stiff(10,4) = -g*jx/le
114 
115  stiff(3,5) = -kwy
116  stiff(5,5) = krr
117  stiff(9,5) = kwy
118  stiff(11,5) = krz
119 
120  stiff(2,6) = kyz
121  stiff(6,6) = kzz
122  stiff(8,6) = -kyz
123  stiff(12,6) = kryz
124 
125  stiff(1,7) = -ea
126  stiff(7,7) = ea
127 
128  stiff(2,8) = -kyy
129  stiff(6,8) = -kyz
130  stiff(8,8) = kyy
131  stiff(12,8) = -kyz
132 
133  stiff(3,9) = -kww
134  stiff(5,9) = kwy
135  stiff(9,9) = kww
136  stiff(11,9) = kwy
137 
138  stiff(4,10) = -g*jx/le
139  stiff(10,10) = g*jx/le
140 
141  stiff(3,11) = -kwy
142  stiff(5,11) = krz
143  stiff(9,11) = kwy
144  stiff(11,11) = krr
145 
146  stiff(2,12) = kyz
147  stiff(6,12) = kryz
148  stiff(8,12) = -kyz
149  stiff(12,12) = kzz
150  end subroutine stf_beam_local
151 
153  subroutine stf_beam(etype,nn,ecoord,section,E,P,STIFF,formulation)
154  integer, intent(in) :: etype
155  integer, intent(in) :: nn
156  real(kind=kreal), intent(in) :: ecoord(3,nn)
157  real(kind=kreal), intent(in) :: section(:)
158  real(kind=kreal), intent(in) :: e,p
159  real(kind=kreal), intent(out) :: stiff(nn*6,nn*6)
160  integer(kind=kint), intent(in), optional :: formulation
161 
162  real(kind=kreal) :: le, trans(3,3), refv(3), transt(3,3)
163 
164  refv = section(1:3)
165  call framtr(refv, ecoord, le, trans)
166  transt= transpose(trans)
167 
168  call stf_beam_local(le, section, e, p, formulation, stiff)
169 
170  stiff(1:3,:) = matmul( transt, stiff(1:3,:) )
171  stiff(4:6,:) = matmul( transt, stiff(4:6,:) )
172  stiff(7:9,:) = matmul( transt, stiff(7:9,:) )
173  stiff(10:12,:) = matmul( transt, stiff(10:12,:) )
174 
175  stiff(:,1:3) = matmul( stiff(:,1:3), trans )
176  stiff(:,4:6) = matmul( stiff(:,4:6), trans )
177  stiff(:,7:9) = matmul( stiff(:,7:9), trans )
178  stiff(:,10:12) = matmul( stiff(:,10:12), trans )
179 
180  end subroutine stf_beam
181 
182  !####################################################################
183  subroutine updatest_beam(etype,nn,ecoord,u,du,section,gausses,QF,formulation)
184  integer, intent(in) :: etype
185  integer, intent(in) :: nn
186  real(kind=kreal), intent(in) :: ecoord(3,nn)
187  real(kind=kreal), intent(in) :: u(6,nn)
188  real(kind=kreal), intent(in) :: du(6,nn)
189  real(kind=kreal), intent(in) :: section(:)
190  type(tgaussstatus), intent(in) :: gausses(:)
191  real(kind=kreal), intent(out) :: qf(nn*6)
192  integer(kind=kint), intent(in), optional :: formulation
193 
194  real(kind=kreal) :: stiff(nn*6, nn*6), totaldisp(nn*6)
195  integer(kind=kint) :: i, j
196  real(kind=kreal) :: e,p
197 
198  e = gausses(1)%pMaterial%variables(m_youngs)
199  p = gausses(1)%pMaterial%variables(m_poisson)
200 
201  call stf_beam(etype,nn,ecoord,section,e,p,stiff,formulation)
202 
203  do i=1,nn
204  do j=1,6
205  totaldisp(6*(i-1)+j) = u(j,i) + du(j,i)
206  end do
207  end do
208 
209  qf = matmul(stiff,totaldisp)
210 
211  end subroutine updatest_beam
212 
213  subroutine stf_beam_641_from_611(ecoord, gausses, section, stiff, formulation)
214  use mmechgauss
215  implicit none
216 
217  real(kind=kreal), intent(in) :: ecoord(3, 4), section(:)
218  type(tgaussstatus), intent(in) :: gausses(:)
219  real(kind=kreal), intent(out) :: stiff(12, 12)
220  integer(kind=kint), intent(in) :: formulation
221 
222  integer(kind=kint) :: i, j
223  integer(kind=kint), parameter :: mixed_to_natural(12) = &
224  (/ 1, 2, 3, 7, 8, 9, 4, 5, 6, 10, 11, 12 /)
225  real(kind=kreal) :: natural_stiff(12, 12), ee, pp
226 
227  ee = gausses(1)%pMaterial%variables(m_youngs)
228  pp = gausses(1)%pMaterial%variables(m_poisson)
229  call stf_beam(611, 2, ecoord(1:3,1:2), section, ee, pp, natural_stiff, formulation)
230 
231  do j = 1, 12
232  do i = 1, 12
233  stiff(i,j) = natural_stiff(mixed_to_natural(i), mixed_to_natural(j))
234  enddo
235  enddo
236  end subroutine stf_beam_641_from_611
237 
238  subroutine updatest_beam_641_from_611(ecoord, u, du, gausses, section, qf, formulation)
239  use mmechgauss
240  implicit none
241 
242  real(kind=kreal), intent(in) :: ecoord(3, 4), u(3, 4), du(3, 4), section(:)
243  type(tgaussstatus), intent(in) :: gausses(:)
244  real(kind=kreal), intent(out) :: qf(12)
245  integer(kind=kint), intent(in) :: formulation
246 
247  integer(kind=kint) :: i
248  integer(kind=kint), parameter :: mixed_to_natural(12) = &
249  (/ 1, 2, 3, 7, 8, 9, 4, 5, 6, 10, 11, 12 /)
250  real(kind=kreal) :: beam_u(6, 2), beam_du(6, 2), natural_qf(12)
251 
252  beam_u(1:3,1:2) = u(1:3,1:2)
253  beam_u(4:6,1:2) = u(1:3,3:4)
254  beam_du(1:3,1:2) = du(1:3,1:2)
255  beam_du(4:6,1:2) = du(1:3,3:4)
256  call updatest_beam(611, 2, ecoord(1:3,1:2), beam_u, beam_du, section, &
257  gausses, natural_qf, formulation)
258 
259  do i = 1, 12
260  qf(i) = natural_qf(mixed_to_natural(i))
261  enddo
262  end subroutine updatest_beam_641_from_611
263 
265  subroutine nqm_beam(nn, ecoord, gausses, section, ul, rnqm, formulation)
266  use mmechgauss
267  integer(kind=kint), intent(in) :: nn
268  real(kind=kreal), intent(in) :: ecoord(3, nn)
269  type(tgaussstatus), intent(in) :: gausses(:)
270  real(kind=kreal), intent(in) :: section(:)
271  real(kind=kreal), intent(in) :: ul(nn*6)
272  real(kind=kreal), intent(out) :: rnqm(nn*6)
273  integer(kind=kint), intent(in), optional :: formulation
274 
275  real(kind=kreal) :: ee, pp, le, refv(3), trans(3,3), ec(3,2)
276  real(kind=kreal) :: stiff(nn*6, nn*6)
277 
278  ee = gausses(1)%pMaterial%variables(m_youngs)
279  pp = gausses(1)%pMaterial%variables(m_poisson)
280 
281  refv(1:3) = section(1:3)
282  ec(1:3, 1) = ecoord(1:3, 1)
283  ec(1:3, 2) = ecoord(1:3, 2)
284  call framtr(refv, ec, le, trans)
285 
286  call stf_beam_local(le, section, ee, pp, formulation, stiff)
287 
288  rnqm = matmul(stiff, ul)
289 
290  end subroutine nqm_beam
291 
293  subroutine nodalstress_beam(etype, nn, ecoord, gausses, section, edisp, ndstrain, ndstress, formulation)
294  use mmechgauss
295  integer(kind=kint), intent(in) :: etype
296  integer(kind=kint), intent(in) :: nn
297  real(kind=kreal), intent(in) :: ecoord(3, nn)
298  type(tgaussstatus), intent(inout) :: gausses(:)
299  real(kind=kreal), intent(in) :: section(:)
300  real(kind=kreal), intent(in) :: edisp(6, nn)
301  real(kind=kreal), intent(out) :: ndstrain(nn, 6)
302  real(kind=kreal), intent(out) :: ndstress(nn, 6)
303  integer(kind=kint), intent(in), optional :: formulation
304 
305  integer(kind=kint) :: k
306  real(kind=kreal) :: ee, pi
307  real(kind=kreal) :: radius, angle(6)
308  real(kind=kreal) :: refv(3), ec(3,2), trans(3,3), le
309  real(kind=kreal) :: ul(nn*6), rnqm(nn*6)
310  real(kind=kreal) :: x2_hat, x3_hat, eps
311  real(kind=kreal) :: stress_i, stress_j
312 
313  pi = 4.0d0*datan(1.0d0)
314 
315  ee = gausses(1)%pMaterial%variables(m_youngs)
316 
317  refv(1:3) = section(1:3)
318  ec(1:3, 1) = ecoord(1:3, 1)
319  ec(1:3, 2) = ecoord(1:3, 2)
320  call framtr(refv, ec, le, trans)
321 
322  radius = gausses(1)%pMaterial%variables(m_beam_radius)
323  angle(1) = gausses(1)%pMaterial%variables(m_beam_angle1)
324  angle(2) = gausses(1)%pMaterial%variables(m_beam_angle2)
325  angle(3) = gausses(1)%pMaterial%variables(m_beam_angle3)
326  angle(4) = gausses(1)%pMaterial%variables(m_beam_angle4)
327  angle(5) = gausses(1)%pMaterial%variables(m_beam_angle5)
328  angle(6) = gausses(1)%pMaterial%variables(m_beam_angle6)
329 
330 
331  ul(1:3) = matmul(trans, edisp(1:3, 1))
332  ul(4:6) = matmul(trans, edisp(4:6, 1))
333  ul(7:9) = matmul(trans, edisp(1:3, 2))
334  ul(10:12) = matmul(trans, edisp(4:6, 2))
335 
336  eps = (ul(7)-ul(1))/le
337  call nqm_beam(nn, ecoord, gausses, section, ul, rnqm, formulation)
338 
339  ndstrain = 0.0d0
340  ndstress = 0.0d0
341 
342  do k = 1, 6
343 
344  angle(k) = angle(k)/180.0d0*pi
345  x2_hat = radius*dcos(angle(k))
346  x3_hat = radius*dsin(angle(k))
347 
348  stress_i = ee*eps + x2_hat*rnqm(6)/section(6) - x3_hat*rnqm(5)/section(5)
349  stress_j = ee*eps - x2_hat*rnqm(12)/section(6) + x3_hat*rnqm(11)/section(5)
350 
351 
352  gausses(1)%strain(k) = eps
353  gausses(1)%stress(k) = 0.5d0*(stress_i + stress_j)
354  gausses(1)%strain_out(k) = gausses(1)%strain(k)
355  gausses(1)%stress_out(k) = gausses(1)%stress(k)
356 
357 
358  ndstrain(1, k) = eps
359  ndstress(1, k) = stress_i
360 
361 
362  ndstrain(2, k) = eps
363  ndstress(2, k) = stress_j
364 
365  end do
366 
367 
368  gausses(1)%nqm(1:nn*6) = rnqm(1:nn*6)
369 
370  end subroutine nodalstress_beam
371 
373  subroutine elementalstress_beam(gausses, estrain, estress, enqm)
374  use mmechgauss
375  type(tgaussstatus), intent(in) :: gausses(:)
376  real(kind=kreal), intent(out) :: estrain(6)
377  real(kind=kreal), intent(out) :: estress(6)
378  real(kind=kreal), intent(out) :: enqm(12)
379 
380  estrain(1:6) = gausses(1)%strain_out(1:6)
381  estress(1:6) = gausses(1)%stress_out(1:6)
382  enqm(1:12) = gausses(1)%nqm(1:12)
383 
384  end subroutine elementalstress_beam
385 
386  ! (Gaku Hashimoto, The University of Tokyo, 2014/02/06) <
388  !####################################################################
389  subroutine stf_beam_641 &
390  (etype, nn, ecoord, gausses, section, stiff, tt, t0)
391  !####################################################################
392 
393  use mmechgauss
394 
395  !--------------------------------------------------------------------
396 
397  integer, intent(in) :: etype
398  integer, intent(in) :: nn
399  real(kind=kreal), intent(in) :: ecoord(3, nn)
400  type(tgaussstatus), intent(in) :: gausses(:)
401  real(kind=kreal), intent(in) :: section(:)
402  real(kind=kreal), intent(out) :: stiff(nn*3, nn*3)
403  real(kind=kreal), intent(in), optional :: tt(nn), t0(nn)
404 
405  !--------------------------------------------------------------------
406 
407  real(kind = kreal) :: refv(3)
408  real(kind = kreal) :: trans(3, 3), transt(3, 3)
409  real(kind = kreal) :: ec(3, 2)
410  real(kind = kreal) :: tempc
411  real(kind = kreal) :: ina1(1), outa1(2)
412  real(kind = kreal) :: ee, pp
413  real(kind = kreal) :: le
414  real(kind = kreal) :: l2, l3, g, a, iy, iz, jx
415  real(kind = kreal) :: ea, twoe, foure, twelvee, sixe
416 
417  logical :: ierr
418 
419  !--------------------------------------------------------------------
420 
421  refv(1) = section(1)
422  refv(2) = section(2)
423  refv(3) = section(3)
424 
425  ec(1, 1) = ecoord(1, 1)
426  ec(2, 1) = ecoord(2, 1)
427  ec(3, 1) = ecoord(3, 1)
428  ec(1, 2) = ecoord(1, 2)
429  ec(2, 2) = ecoord(2, 2)
430  ec(3, 2) = ecoord(3, 2)
431 
432  call framtr(refv, ec, le, trans)
433 
434  transt= transpose( trans )
435 
436  l2 = le*le
437  l3 = l2*le
438 
439  !--------------------------------------------------------------------
440 
441  if( present( tt ) ) then
442 
443  tempc = 0.5d0*( tt(1)+tt(2) )
444 
445  end if
446 
447  !--------------------------------------------------------------------
448 
449  if( present( tt ) ) then
450 
451  ina1(1) = tempc
452 
453  call fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa1, ierr, ina1 )
454 
455  else
456 
457  ierr = .true.
458 
459  end if
460 
461  !--------------------------------------------------------------
462 
463  if( ierr ) then
464 
465  ee = gausses(1)%pMaterial%variables(m_youngs)
466  pp = gausses(1)%pMaterial%variables(m_poisson)
467 
468  else
469 
470  ee = outa1(1)
471  pp = outa1(2)
472 
473  end if
474 
475  !--------------------------------------------------------------------
476 
477  g = ee/( 2.0d0*( 1.0d0+pp ) )
478 
479  a = section(4)
480 
481  iy = section(5)
482  iz = section(6)
483  jx = section(7)
484 
485  !--------------------------------------------------------------------
486 
487  ea = ee*a/le
488 
489  twoe = 2.0d0*ee/le
490  foure = 4.0d0*ee/le
491  twelvee = 12.0d0*ee/l3
492  sixe = 6.0d0*ee/l2
493 
494  !--------------------------------------------------------------------
495 
496  stiff = 0.0d0
497 
498  stiff(1, 1) = ea
499  !stiff(7, 1) = -ea
500  stiff(4, 1) = -ea
501 
502  stiff(2, 2) = twelvee*iz
503  !stiff(6, 2) = sixe*iz
504  stiff(9, 2) = sixe*iz
505  !stiff(8, 2) = -twelvee*iz
506  stiff(5, 2) = -twelvee*iz
507  stiff(12, 2) = sixe*iz
508 
509  stiff(3, 3) = twelvee*iy
510  !stiff(5, 3) = -sixe*iy
511  stiff(8, 3) = -sixe*iy
512  !stiff(9, 3) = -twelvee*iy
513  stiff(6, 3) = -twelvee*iy
514  stiff(11, 3) = -sixe*iy
515 
516  !stiff(4, 4) = g*jx/le
517  stiff(7, 7) = g*jx/le
518  !stiff(10, 4) = -g*jx/le
519  stiff(10, 7) = -g*jx/le
520 
521  !stiff(3, 5) = -sixe*iy
522  stiff(3, 8) = -sixe*iy
523  !stiff(5, 5) = foure*iy
524  stiff(8, 8) = foure*iy
525  !stiff(9, 5) = sixe*iy
526  stiff(6, 8) = sixe*iy
527  !stiff(11, 5) = twoe*iy
528  stiff(11, 8) = twoe*iy
529 
530  !stiff(2, 6) = sixe*iz
531  stiff(2, 9) = sixe*iz
532  !stiff(6, 6) = foure*iz
533  stiff(9, 9) = foure*iz
534  !stiff(8, 6) = -sixe*iz
535  stiff(5, 9) = -sixe*iz
536  !stiff(12, 6) = twoe*iz
537  stiff(12, 9) = twoe*iz
538 
539  !stiff(1, 7) = -ea
540  stiff(1, 4) = -ea
541  !stiff(7, 7) = ea
542  stiff(4, 4) = ea
543 
544  !stiff(2, 8) = -twelvee*iz
545  stiff(2, 5) = -twelvee*iz
546  !stiff(6, 8) = -sixe*iz
547  stiff(9, 5) = -sixe*iz
548  !stiff(8, 8) = twelvee*iz
549  stiff(5, 5) = twelvee*iz
550  !stiff(12, 8) = -sixe*iz
551  stiff(12, 5) = -sixe*iz
552 
553  !stiff(3, 9) = -twelvee*iy
554  stiff(3, 6) = -twelvee*iy
555  !stiff(5, 9) = sixe*iy
556  stiff(8, 6) = sixe*iy
557  !stiff(9, 9) = twelvee*iy
558  stiff(6, 6) = twelvee*iy
559  !stiff(11, 9) = sixe*iy
560  stiff(11, 6) = sixe*iy
561 
562  !stiff(4, 10) = -g*jx/le
563  stiff(7, 10) = -g*jx/le
564  stiff(10, 10) = g*jx/le
565 
566  stiff(3, 11) = -sixe*iy
567  !stiff(5, 11) = twoe*iy
568  stiff(8, 11) = twoe*iy
569  !stiff(9, 11) = sixe*iy
570  stiff(6, 11) = sixe*iy
571  stiff(11, 11) = foure*iy
572 
573  stiff(2, 12) = sixe*iz
574  !stiff(6, 12) = twoe*iz
575  stiff(9, 12) = twoe*iz
576  !stiff(8, 12) = -sixe*iz
577  stiff(5, 12) = -sixe*iz
578  stiff(12, 12) = foure*iz
579 
580  !--------------------------------------------------------------------
581 
582  stiff( 1:3, :) = matmul( transt, stiff( 1:3, :) )
583  stiff( 4:6, :) = matmul( transt, stiff( 4:6, :) )
584  stiff( 7:9, :) = matmul( transt, stiff( 7:9, :) )
585  stiff(10:12, :) = matmul( transt, stiff(10:12, :) )
586 
587  stiff(:, 1:3) = matmul( stiff(:, 1:3), trans )
588  stiff(:, 4:6) = matmul( stiff(:, 4:6), trans )
589  stiff(:, 7:9) = matmul( stiff(:, 7:9), trans )
590  stiff(:, 10:12) = matmul( stiff(:, 10:12), trans )
591 
592  !--------------------------------------------------------------------
593 
594  return
595 
596  !####################################################################
597  end subroutine stf_beam_641
598  !####################################################################
599  ! > (Gaku Hashimoto, The University of Tokyo, 2014/02/06)
600 
602 !####################################################################
603  SUBROUTINE nqm_beam_641 &
604  (etype, nn, ecoord, gausses, section, stiff, tt, t0, tdisp, rnqm)
605 !####################################################################
606 
607  USE mmechgauss
608 
609 !--------------------------------------------------------------------
610 
611  INTEGER, INTENT(IN) :: etype
612  INTEGER, INTENT(IN) :: nn
613  REAL(kind=kreal), INTENT(IN) :: ecoord(3, nn)
614  TYPE(tgaussstatus), INTENT(IN) :: gausses(:)
615  REAL(kind=kreal), INTENT(IN) :: section(:)
616  REAL(kind=kreal), INTENT(OUT) :: stiff(nn*3, nn*3)
617  REAL(kind=kreal), INTENT(IN), OPTIONAL :: tt(nn), t0(nn)
618 
619  REAL(kind=kreal), INTENT(INOUT) :: tdisp(nn*3)
620  REAL(kind=kreal), INTENT(OUT) :: rnqm(nn*3)
621 
622  REAL(kind=kreal) :: tdisp1(nn*3)
623 
624 !--------------------------------------------------------------------
625 
626  REAL(kind = kreal) :: refv(3)
627  REAL(kind = kreal) :: trans(3, 3), transt(3, 3)
628  REAL(kind = kreal) :: ec(3, 2)
629  REAL(kind = kreal) :: tempc
630  REAL(kind = kreal) :: ina1(1), outa1(2)
631  REAL(kind = kreal) :: ee, pp
632  REAL(kind = kreal) :: le
633  REAL(kind = kreal) :: l2, l3, g, a, iy, iz, jx
634  REAL(kind = kreal) :: ea, twoe, foure, twelvee, sixe
635 
636  LOGICAL :: ierr
637 
638 !--------------------------------------------------------------------
639 
640  refv(1) = section(1)
641  refv(2) = section(2)
642  refv(3) = section(3)
643 
644  ec(1, 1) = ecoord(1, 1)
645  ec(2, 1) = ecoord(2, 1)
646  ec(3, 1) = ecoord(3, 1)
647  ec(1, 2) = ecoord(1, 2)
648  ec(2, 2) = ecoord(2, 2)
649  ec(3, 2) = ecoord(3, 2)
650 
651  CALL framtr(refv, ec, le, trans)
652 
653  transt= transpose( trans )
654 
655  l2 = le*le
656  l3 = l2*le
657 
658 !--------------------------------------------------------------------
659 
660  IF( PRESENT( tt ) ) THEN
661 
662  tempc = 0.5d0*( tt(1)+tt(2) )
663 
664  END IF
665 
666 !--------------------------------------------------------------------
667 
668  IF( PRESENT( tt ) ) THEN
669 
670  ina1(1) = tempc
671 
672  CALL fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa1, ierr, ina1 )
673 
674  ELSE
675 
676  ierr = .true.
677 
678  END IF
679 
680  !--------------------------------------------------------------
681 
682  IF( ierr ) THEN
683 
684  ee = gausses(1)%pMaterial%variables(m_youngs)
685  pp = gausses(1)%pMaterial%variables(m_poisson)
686 
687  ELSE
688 
689  ee = outa1(1)
690  pp = outa1(2)
691 
692  END IF
693 
694 !--------------------------------------------------------------------
695 
696  g = ee/( 2.0d0*( 1.0d0+pp ) )
697 
698  a = section(4)
699 
700  iy = section(5)
701  iz = section(6)
702  jx = section(7)
703 
704 ! write (6,'(a,4e15.5)') 'a,iy,iz,jx',a,iy,iz,jx
705 
706 !--------------------------------------------------------------------
707 
708  ea = ee*a/le
709 
710  twoe = 2.0d0*ee/le
711  foure = 4.0d0*ee/le
712  twelvee = 12.0d0*ee/l3
713  sixe = 6.0d0*ee/l2
714 
715 !--------------------------------------------------------------------
716 
717  stiff = 0.0d0
718 
719  stiff(1, 1) = ea
720  !stiff(7, 1) = -ea
721  stiff(4, 1) = -ea
722 
723  stiff(2, 2) = twelvee*iz
724  !stiff(6, 2) = sixe*iz
725  stiff(9, 2) = sixe*iz
726  !stiff(8, 2) = -twelvee*iz
727  stiff(5, 2) = -twelvee*iz
728  stiff(12, 2) = sixe*iz
729 
730  stiff(3, 3) = twelvee*iy
731  !stiff(5, 3) = -sixe*iy
732  stiff(8, 3) = -sixe*iy
733  !stiff(9, 3) = -twelvee*iy
734  stiff(6, 3) = -twelvee*iy
735  stiff(11, 3) = -sixe*iy
736 
737  !stiff(4, 4) = g*jx/le
738  stiff(7, 7) = g*jx/le
739  !stiff(10, 4) = -g*jx/le
740  stiff(10, 7) = -g*jx/le
741 
742  !stiff(3, 5) = -sixe*iy
743  stiff(3, 8) = -sixe*iy
744  !stiff(5, 5) = foure*iy
745  stiff(8, 8) = foure*iy
746  !stiff(9, 5) = sixe*iy
747  stiff(6, 8) = sixe*iy
748  !stiff(11, 5) = twoe*iy
749  stiff(11, 8) = twoe*iy
750 
751  !stiff(2, 6) = sixe*iz
752  stiff(2, 9) = sixe*iz
753  !stiff(6, 6) = foure*iz
754  stiff(9, 9) = foure*iz
755  !stiff(8, 6) = -sixe*iz
756  stiff(5, 9) = -sixe*iz
757  !stiff(12, 6) = twoe*iz
758  stiff(12, 9) = twoe*iz
759 
760  !stiff(1, 7) = -ea
761  stiff(1, 4) = -ea
762  !stiff(7, 7) = ea
763  stiff(4, 4) = ea
764 
765  !stiff(2, 8) = -twelvee*iz
766  stiff(2, 5) = -twelvee*iz
767  !stiff(6, 8) = -sixe*iz
768  stiff(9, 5) = -sixe*iz
769  !stiff(8, 8) = twelvee*iz
770  stiff(5, 5) = twelvee*iz
771  !stiff(12, 8) = -sixe*iz
772  stiff(12, 5) = -sixe*iz
773 
774  !stiff(3, 9) = -twelvee*iy
775  stiff(3, 6) = -twelvee*iy
776  !stiff(5, 9) = sixe*iy
777  stiff(8, 6) = sixe*iy
778  !stiff(9, 9) = twelvee*iy
779  stiff(6, 6) = twelvee*iy
780  !stiff(11, 9) = sixe*iy
781  stiff(11, 6) = sixe*iy
782 
783  !stiff(4, 10) = -g*jx/le
784  stiff(7, 10) = -g*jx/le
785  stiff(10, 10) = g*jx/le
786 
787  stiff(3, 11) = -sixe*iy
788  !stiff(5, 11) = twoe*iy
789  stiff(8, 11) = twoe*iy
790  !stiff(9, 11) = sixe*iy
791  stiff(6, 11) = sixe*iy
792  stiff(11, 11) = foure*iy
793 
794  stiff(2, 12) = sixe*iz
795  !stiff(6, 12) = twoe*iz
796  stiff(9, 12) = twoe*iz
797  !stiff(8, 12) = -sixe*iz
798  stiff(5, 12) = -sixe*iz
799  stiff(12, 12) = foure*iz
800 
801 !--------------------------------------------------------------------
802  tdisp1( 1: 3 ) = matmul( trans, tdisp( 1: 3 ) )
803  tdisp1( 4: 6 ) = matmul( trans, tdisp( 4: 6 ) )
804  tdisp1( 7: 9 ) = matmul( trans, tdisp( 7: 9 ) )
805  tdisp1( 10:12 ) = matmul( trans, tdisp( 10:12 ) )
806 !--------------------------------------------------------------------
807  rnqm( 1:12 ) = matmul( stiff, tdisp1 )
808 
809  RETURN
810 
811 !####################################################################
812  END SUBROUTINE nqm_beam_641
813 !####################################################################
814 
815  ! (Gaku Hashimoto, The University of Tokyo, 2014/02/06) <
816  !####################################################################
817  subroutine dl_beam_641(etype, nn, xx, yy, zz, rho, ltype, params, &
818  section, vect, nsize)
819  !####################################################################
820  !**
821  !** SET DLOAD
822  !**
823  ! BX LTYPE=1 :BODY FORCE IN X-DIRECTION
824  ! BY LTYPE=2 :BODY FORCE IN Y-DIRECTION
825  ! BZ LTYPE=3 :BODY FORCE IN Z-DIRECTION
826  ! GRAV LTYPE=4 :GRAVITY FORCE
827  ! CENT LTYPE=5 :CENTRIFUGAL LOAD
828  ! P1 LTYPE=10 :TRACTION IN NORMAL-DIRECTION FOR FACE-1
829  ! P2 LTYPE=20 :TRACTION IN NORMAL-DIRECTION FOR FACE-2
830  ! P3 LTYPE=30 :TRACTION IN NORMAL-DIRECTION FOR FACE-3
831  ! P4 LTYPE=40 :TRACTION IN NORMAL-DIRECTION FOR FACE-4
832  ! P5 LTYPE=50 :TRACTION IN NORMAL-DIRECTION FOR FACE-5
833  ! P6 LTYPE=60 :TRACTION IN NORMAL-DIRECTION FOR FACE-6
834  ! I/F VARIABLES
835  integer(kind = kint), intent(in) :: etype, nn
836  real(kind = kreal), intent(in) :: xx(:), yy(:), zz(:)
837  real(kind = kreal), intent(in) :: params(0:6)
838  real(kind = kreal), intent(in) :: section(:)
839  real(kind = kreal), intent(inout) :: vect(:)
840  real(kind = kreal) :: rho
841  integer(kind = kint) :: ltype, nsize
842  ! LOCAL VARIABLES
843  integer(kind = kint) :: ndof
844  parameter(ndof = 3)
845  integer(kind = kint) :: ivol, isuf, nod(nn)
846  integer(kind = kint) :: i ,surtype, nsur
847  real(kind = kreal) :: vx, vy, vz, val, a, aa
848 
849  !--------------------------------------------------------------------
850 
851  val = params(0)
852 
853  !--------------------------------------------------------------
854 
855  ivol = 0
856  isuf = 0
857 
858  if( ltype .LT. 10 ) then
859 
860  ivol = 1
861 
862  else if( ltype .GE. 10 ) then
863 
864  isuf = 1
865 
866  call getsubface(etype, ltype/10, surtype, nod)
867 
868  nsur = getnumberofnodes(surtype)
869 
870  end if
871 
872  !--------------------------------------------------------------------
873 
874  nsize = nn*ndof
875 
876  !--------------------------------------------------------------------
877 
878  vect(1:nsize) = 0.0d0
879 
880  !--------------------------------------------------------------
881 
882  ! Volume force
883 
884  if( ivol .EQ. 1 ) then
885 
886  if( ltype .EQ. 4 ) then
887 
888  aa = dsqrt( ( xx(2)-xx(1) )*( xx(2)-xx(1) ) &
889  +( yy(2)-yy(1) )*( yy(2)-yy(1) ) &
890  +( zz(2)-zz(1) )*( zz(2)-zz(1) ) )
891 
892  a = section(4)
893 
894  vx = params(1)
895  vy = params(2)
896  vz = params(3)
897  vx = vx/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
898  vy = vy/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
899  vz = vz/dsqrt( params(1)**2+params(2)**2+params(3)**2 )
900 
901  do i = 1, 2
902 
903  vect(3*i-2) = val*rho*a*0.5d0*aa*vx
904  vect(3*i-1) = val*rho*a*0.5d0*aa*vy
905  vect(3*i ) = val*rho*a*0.5d0*aa*vz
906 
907  end do
908 
909  do i = 3, 4
910 
911  vect(3*i-2) = 0.0d0
912  vect(3*i-1) = 0.0d0
913  vect(3*i ) = 0.0d0
914 
915  end do
916 
917  end if
918 
919  end if
920 
921  !--------------------------------------------------------------------
922 
923  return
924 
925  !####################################################################
926  end subroutine dl_beam_641
927  !####################################################################
928  ! > (Gaku Hashimoto, The University of Tokyo, 2014/02/06)
929 
930 
931  ! (Gaku Hashimoto, The University of Tokyo, 2014/02/06) <
932  !####################################################################
933  subroutine tload_beam_641 &
934  (etype, nn, ndof, xx, yy, zz, tt, t0, &
935  gausses, section, vect)
936  !####################################################################
937 
938  use hecmw
939  use m_fstr
940  use m_utilities
941  use mmechgauss
942 
943  !--------------------------------------------------------------------
944 
945  integer(kind = kint), intent(in) :: etype
946  integer(kind = kint), intent(in) :: nn
947  integer(kind = kint), intent(in) :: ndof
948  type(tgaussstatus), intent(in) :: gausses(:)
949  real(kind = kreal), intent(in) :: section(:)
950  real(kind = kreal), intent(in) :: xx(nn), yy(nn), zz(nn)
951  real(kind = kreal), intent(in) :: tt(nn), t0(nn)
952  real(kind = kreal), intent(out) :: vect(nn*ndof)
953 
954  !--------------------------------------------------------------------
955 
956  real(kind = kreal) :: tempc, temp0
957  real(kind = kreal) :: ecoord(3, nn)
958  real(kind = kreal) :: ec(3, 2)
959  real(kind = kreal) :: ina1(1), outa1(2)
960  real(kind = kreal) :: ina2(1), outa2(1)
961  real(kind = kreal) :: alp, alp0
962  real(kind = kreal) :: ee, pp
963  real(kind = kreal) :: a
964  real(kind = kreal) :: refv(3)
965  real(kind = kreal) :: g
966  real(kind = kreal) :: le
967  real(kind = kreal) :: trans(3, 3), transt(3, 3)
968 
969  logical :: ierr
970 
971  !--------------------------------------------------------------------
972 
973  ecoord(1, 1:nn) = xx(1:nn)
974  ecoord(2, 1:nn) = yy(1:nn)
975  ecoord(3, 1:nn) = zz(1:nn)
976 
977  !--------------------------------------------------------------------
978 
979  tempc = 0.5d0*( tt(1)+tt(2) )
980  temp0 = 0.5d0*( t0(1)+t0(2) )
981 
982  !--------------------------------------------------------------
983 
984  ina1(1) = tempc
985 
986  call fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa1, ierr, ina1 )
987 
988  if( ierr ) then
989 
990  ee = gausses(1)%pMaterial%variables(m_youngs)
991  pp = gausses(1)%pMaterial%variables(m_poisson)
992 
993  else
994 
995  ee = outa1(1)
996  pp = outa1(2)
997 
998  end if
999 
1000  !--------------------------------------------------------------
1001 
1002  ina2(1) = tempc
1003 
1004  call fetch_tabledata( mc_themoexp, gausses(1)%pMaterial%dict, outa2(:), ierr, ina2 )
1005 
1006  if( ierr ) stop "Fails in fetching expansion coefficient!"
1007 
1008  alp = outa2(1)
1009 
1010  !--------------------------------------------------------------
1011 
1012  ina2(1) = temp0
1013 
1014  call fetch_tabledata( mc_themoexp, gausses(1)%pMaterial%dict, outa2(:), ierr, ina2 )
1015 
1016  if( ierr ) stop "Fails in fetching expansion coefficient!"
1017 
1018  alp0 = outa2(1)
1019 
1020  !--------------------------------------------------------------------
1021 
1022  refv(1) = section(1)
1023  refv(2) = section(2)
1024  refv(3) = section(3)
1025 
1026  ec(1, 1) = ecoord(1, 1)
1027  ec(2, 1) = ecoord(2, 1)
1028  ec(3, 1) = ecoord(3, 1)
1029  ec(1, 2) = ecoord(1, 2)
1030  ec(2, 2) = ecoord(2, 2)
1031  ec(3, 2) = ecoord(3, 2)
1032 
1033  call framtr(refv, ec, le, trans)
1034 
1035  transt= transpose( trans )
1036 
1037  !--------------------------------------------------------------------
1038 
1039  a = section(4)
1040 
1041  g = ee/( 2.0d0*( 1.0d0+pp ))
1042 
1043  !--------------------------------------------------------------------
1044 
1045  vect( 1) = -a*ee*( alp*( tempc-ref_temp )-alp0*( temp0-ref_temp ) )
1046  vect( 2) = 0.0d0
1047  vect( 3) = 0.0d0
1048 
1049  vect( 4) = a*ee*( alp*( tempc-ref_temp )-alp0*( temp0-ref_temp ) )
1050  vect( 5) = 0.0d0
1051  vect( 6) = 0.0d0
1052 
1053  vect( 7) = 0.0d0
1054  vect( 8) = 0.0d0
1055  vect( 9) = 0.0d0
1056 
1057  vect(10) = 0.0d0
1058  vect(11) = 0.0d0
1059  vect(12) = 0.0d0
1060 
1061  !--------------------------------------------------------------------
1062 
1063  vect( 1:3) = matmul( transt, vect(1:3) )
1064  vect( 4:6) = matmul( transt, vect(4:6) )
1065  vect( 7:9) = matmul( transt, vect(7:9) )
1066  vect(10:12) = matmul( transt, vect(10:12) )
1067 
1068  !--------------------------------------------------------------------
1069 
1070  return
1071 
1072  !####################################################################
1073  end subroutine tload_beam_641
1074  !####################################################################
1075  ! > (Gaku Hashimoto, The University of Tokyo, 2013/09/13)
1076 
1077 
1078  ! (Gaku Hashimoto, The University of Tokyo, 2014/02/06) <
1079  !####################################################################
1080  subroutine nodalstress_beam_641 &
1081  (etype, nn, ecoord, gausses, section, edisp, &
1082  ndstrain, ndstress, tt, t0, ntemp)
1083  !####################################################################
1084 
1085  use m_fstr
1086  use mmechgauss
1087 
1088  !--------------------------------------------------------------------
1089 
1090  integer(kind = kint), intent(in) :: etype
1091  integer(kind = kint), intent(in) :: nn
1092  real(kind = kreal), intent(in) :: ecoord(3, nn)
1093  type(tgaussstatus), intent(inout) :: gausses(:)
1094  real(kind = kreal), intent(in) :: section(:)
1095  real(kind = kreal), intent(in) :: edisp(3, nn)
1096  real(kind = kreal), intent(out) :: ndstrain(nn, 6)
1097  real(kind = kreal), intent(out) :: ndstress(nn, 6)
1098  real(kind=kreal), intent(in), optional :: tt(nn), t0(nn)
1099  integer(kind = kint), intent(in) :: ntemp
1100 
1101  !--------------------------------------------------------------------
1102 
1103  real(kind=kreal) :: stiffx(12, 12)
1104  real(kind=kreal) :: tdisp(12)
1105  real(kind=kreal) :: rnqm(12)
1106 
1107  !--------------------------------------------------------------------
1108 
1109  integer(kind = kint) :: i, j, k, jj
1110 
1111  real(kind = kreal) :: tempc, temp0
1112  real(kind = kreal) :: ina1(1), outa1(2)
1113  real(kind = kreal) :: ina2(1), outa2(1)
1114  real(kind = kreal) :: alp, alp0
1115  real(kind = kreal) :: ee, pp
1116  real(kind = kreal) :: a, radius, angle(6)
1117  real(kind = kreal) :: refv(3)
1118  real(kind = kreal) :: le, l2, l3
1119  real(kind = kreal) :: trans(3, 3), transt(3, 3)
1120  real(kind = kreal) :: edisp_hat(3, nn)
1121  real(kind = kreal) :: ec(3, 2)
1122  real(kind = kreal) :: t(3, 3), t_hat(3, 3)
1123  real(kind = kreal) :: t_hat_tmp(3, 3)
1124  real(kind = kreal) :: e(3, 3), e_hat(3, 3)
1125  real(kind = kreal) :: e_hat_tmp(3, 3)
1126  real(kind = kreal) :: x1_hat, x2_hat, x3_hat
1127  real(kind = kreal) :: pi
1128 
1129  logical :: ierr
1130 
1131  alp = 0.0d0; alp0 = 0.0d0
1132  tempc = 0.0d0; temp0 = 0.0d0
1133 
1134  !--------------------------------------------------------------------
1135 
1136  pi = 4.0d0*datan( 1.0d0 )
1137 
1138  !--------------------------------------------------------------------
1139 
1140  if( present( tt ) .AND. present( t0 ) ) then
1141 
1142  tempc = 0.5d0*( tt(1)+tt(2) )
1143  temp0 = 0.5d0*( t0(1)+t0(2) )
1144 
1145  end if
1146 
1147  !--------------------------------------------------------------------
1148 
1149  if( ntemp .EQ. 1 ) then
1150 
1151  ina1(1) = tempc
1152 
1153  call fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa1, ierr, ina1 )
1154 
1155  else
1156 
1157  ierr = .true.
1158 
1159  end if
1160 
1161  !--------------------------------------------------------------
1162 
1163  if( ierr ) then
1164 
1165  ee = gausses(1)%pMaterial%variables(m_youngs)
1166  pp = gausses(1)%pMaterial%variables(m_poisson)
1167 
1168  else
1169 
1170  ee = outa1(1)
1171  pp = outa1(2)
1172 
1173  end if
1174 
1175  !--------------------------------------------------------------------
1176 
1177  if( ntemp .EQ. 1 ) then
1178 
1179  ina2(1) = tempc
1180 
1181  call fetch_tabledata( mc_themoexp, gausses(1)%pMaterial%dict, outa2(:), ierr, ina2 )
1182 
1183  if( ierr ) stop "Fails in fetching expansion coefficient!"
1184 
1185  alp = outa2(1)
1186 
1187  end if
1188 
1189  !--------------------------------------------------------------
1190 
1191  if( ntemp .EQ. 1 ) then
1192 
1193  ina2(1) = temp0
1194 
1195  call fetch_tabledata( mc_themoexp, gausses(1)%pMaterial%dict, outa2(:), ierr, ina2 )
1196 
1197  if( ierr ) stop "Fails in fetching expansion coefficient!"
1198 
1199  alp0 = outa2(1)
1200 
1201  end if
1202 
1203  !--------------------------------------------------------------------
1204 
1205  refv(1) = section(1)
1206  refv(2) = section(2)
1207  refv(3) = section(3)
1208 
1209  ec(1, 1) = ecoord(1, 1)
1210  ec(2, 1) = ecoord(2, 1)
1211  ec(3, 1) = ecoord(3, 1)
1212  ec(1, 2) = ecoord(1, 2)
1213  ec(2, 2) = ecoord(2, 2)
1214  ec(3, 2) = ecoord(3, 2)
1215 
1216  call framtr(refv, ec, le, trans)
1217 
1218  transt= transpose( trans )
1219 
1220  l2 = le*le
1221  l3 = l2*le
1222 
1223  !--------------------------------------------------------------------
1224 
1225  a = section(4)
1226 
1227  radius = gausses(1)%pMaterial%variables(m_beam_radius)
1228 
1229  angle(1) = gausses(1)%pMaterial%variables(m_beam_angle1)
1230  angle(2) = gausses(1)%pMaterial%variables(m_beam_angle2)
1231  angle(3) = gausses(1)%pMaterial%variables(m_beam_angle3)
1232  angle(4) = gausses(1)%pMaterial%variables(m_beam_angle4)
1233  angle(5) = gausses(1)%pMaterial%variables(m_beam_angle5)
1234  angle(6) = gausses(1)%pMaterial%variables(m_beam_angle6)
1235 
1236  !--------------------------------------------------------------------
1237 
1238  do k = 1, 6
1239 
1240  !--------------------------------------------------------
1241 
1242  angle(k) = angle(k)/180.0d0*pi
1243 
1244  x2_hat = radius*dcos( angle(k) )
1245  x3_hat = radius*dsin( angle(k) )
1246 
1247  !--------------------------------------------------------
1248 
1249  jj = 0
1250  do j = 1, nn
1251 
1252  do i = 1, 3
1253 
1254  edisp_hat(i, j) = trans(i, 1)*edisp(1, j) &
1255  +trans(i, 2)*edisp(2, j) &
1256  +trans(i, 3)*edisp(3, j)
1257 
1258  jj = jj + 1
1259  tdisp(jj) = edisp(i,j)
1260 
1261  end do
1262 
1263  end do
1264 
1265  !--------------------------------------------------------
1266 
1267  x1_hat = 0.5d0*le
1268 
1269  e_hat = 0.0d0
1270  e_hat(1, 1) = ( edisp_hat(1, 2)-edisp_hat(1, 1) )/le
1271 
1272  t_hat = 0.0d0
1273  t_hat(1, 1) = ee*( edisp_hat(1, 2)-edisp_hat(1, 1) )/le &
1274  -ee*x2_hat*( ( -6.0d0/l2+12.0d0*x1_hat/l3 )*edisp_hat(2, 1) &
1275  +( -4.0d0/le+6.0d0*x1_hat/l2 )*edisp_hat(3, 3) &
1276  +( 6.0d0/l2-12.0d0*x1_hat/l3 )*edisp_hat(2, 2) &
1277  +( -2.0d0/le+6.0d0*x1_hat/l2 )*edisp_hat(3, 4) ) &
1278  -ee*x3_hat*( ( -6.0d0/l2+12.0d0*x1_hat/l3 )*edisp_hat(3, 1) &
1279  +( 4.0d0/le-6.0d0*x1_hat/l2 )*edisp_hat(2, 3) &
1280  +( 6.0d0/l2-12.0d0*x1_hat/l3 )*edisp_hat(3, 2) &
1281  +( 2.0d0/le-6.0d0*x1_hat/l2 )*edisp_hat(2, 4) )
1282 
1283  if( ntemp .EQ. 1 ) then
1284 
1285  t_hat(1, 1) &
1286  = t_hat(1, 1) &
1287  -ee*( alp*( tempc-ref_temp )-alp0*( temp0-ref_temp ) )
1288 
1289  end if
1290 
1291  e_hat_tmp(1:3,:) = matmul( trans, e_hat(1:3,:) )
1292  t_hat_tmp(1:3,:) = matmul( trans, t_hat(1:3,:) )
1293 
1294  e(:, 1:3) = matmul( e_hat_tmp(:,1:3), transt )
1295  t(:, 1:3) = matmul( t_hat_tmp(:,1:3), transt )
1296 
1297  gausses(1)%strain(k) = e_hat(1, 1)
1298  gausses(1)%stress(k) = t_hat(1, 1)
1299 
1300  !set stress and strain for output
1301  gausses(1)%strain_out(k) = gausses(1)%strain(k)
1302  gausses(1)%stress_out(k) = gausses(1)%stress(k)
1303 
1304  !--------------------------------------------------------
1305 
1306  ndstrain(1, k) = 0.0d0
1307  ndstrain(2, k) = 0.0d0
1308  ndstrain(3, k) = 0.0d0
1309  ndstrain(4, k) = 0.0d0
1310 
1311  ndstress(1, k) = 0.0d0
1312  ndstress(2, k) = 0.0d0
1313  ndstress(3, k) = 0.0d0
1314  ndstress(4, k) = 0.0d0
1315 
1316  !--------------------------------------------------------
1317 
1318  x1_hat = 0.0d0
1319 
1320  e_hat = 0.0d0
1321  e_hat(1, 1) = ( edisp_hat(1, 2)-edisp_hat(1, 1) )/le
1322 
1323  t_hat = 0.0d0
1324  t_hat(1, 1) = ee*( edisp_hat(1, 2)-edisp_hat(1, 1) )/le &
1325  -ee*x2_hat*( ( -6.0d0/l2+12.0d0*x1_hat/l3 )*edisp_hat(2, 1) &
1326  +( -4.0d0/le+6.0d0*x1_hat/l2 )*edisp_hat(3, 3) &
1327  +( 6.0d0/l2-12.0d0*x1_hat/l3 )*edisp_hat(2, 2) &
1328  +( -2.0d0/le+6.0d0*x1_hat/l2 )*edisp_hat(3, 4) ) &
1329  -ee*x3_hat*( ( -6.0d0/l2+12.0d0*x1_hat/l3 )*edisp_hat(3, 1) &
1330  +( 4.0d0/le-6.0d0*x1_hat/l2 )*edisp_hat(2, 3) &
1331  +( 6.0d0/l2-12.0d0*x1_hat/l3 )*edisp_hat(3, 2) &
1332  +( 2.0d0/le-6.0d0*x1_hat/l2 )*edisp_hat(2, 4) )
1333 
1334  if( ntemp .EQ. 1 ) then
1335 
1336  t_hat(1, 1) &
1337  = t_hat(1, 1) &
1338  -ee*( alp*( tempc-ref_temp )-alp0*( temp0-ref_temp ) )
1339 
1340  end if
1341 
1342  e_hat_tmp(1:3, :) = matmul( trans, e_hat(1:3, :) )
1343  t_hat_tmp(1:3, :) = matmul( trans, t_hat(1:3, :) )
1344 
1345  e(:, 1:3) = matmul( e_hat_tmp(:, 1:3), transt )
1346  t(:, 1:3) = matmul( t_hat_tmp(:, 1:3), transt )
1347 
1348  ndstrain(1, k) = e_hat(1, 1)
1349  ndstress(1, k) = t_hat(1, 1)
1350 
1351  !--------------------------------------------------------
1352 
1353  x1_hat = le
1354 
1355  e_hat = 0.0d0
1356  e_hat(1, 1) = ( edisp_hat(1, 2)-edisp_hat(1, 1) )/le
1357 
1358  t_hat = 0.0d0
1359  t_hat(1, 1) = ee*( edisp_hat(1, 2)-edisp_hat(1, 1) )/le &
1360  -ee*x2_hat*( ( -6.0d0/l2+12.0d0*x1_hat/l3 )*edisp_hat(2, 1) &
1361  +( -4.0d0/le+6.0d0*x1_hat/l2 )*edisp_hat(3, 3) &
1362  +( 6.0d0/l2-12.0d0*x1_hat/l3 )*edisp_hat(2, 2) &
1363  +( -2.0d0/le+6.0d0*x1_hat/l2 )*edisp_hat(3, 4) ) &
1364  -ee*x3_hat*( ( -6.0d0/l2+12.0d0*x1_hat/l3 )*edisp_hat(3, 1) &
1365  +( 4.0d0/le-6.0d0*x1_hat/l2 )*edisp_hat(2, 3) &
1366  +( 6.0d0/l2-12.0d0*x1_hat/l3 )*edisp_hat(3, 2) &
1367  +( 2.0d0/le-6.0d0*x1_hat/l2 )*edisp_hat(2, 4) )
1368 
1369  if( ntemp .EQ. 1 ) then
1370 
1371  t_hat(1, 1) &
1372  = t_hat(1, 1) &
1373  -ee*( alp*( tempc-ref_temp )-alp0*( temp0-ref_temp ) )
1374 
1375  end if
1376 
1377  e_hat_tmp(1:3, :) = matmul( trans, e_hat(1:3, :) )
1378  t_hat_tmp(1:3, :) = matmul( trans, t_hat(1:3, :) )
1379 
1380  e(:, 1:3) = matmul( e_hat_tmp(:, 1:3), transt )
1381  t(:, 1:3) = matmul( t_hat_tmp(:, 1:3), transt )
1382 
1383  ndstrain(2, k) = e_hat(1, 1)
1384  ndstress(2, k) = t_hat(1, 1)
1385 
1386  !--------------------------------------------------------
1387 
1388  end do
1389 
1390  !--------------------------------------------------------------------
1391  stiffx = 0.0
1392 
1393  call nqm_beam_641 &
1394  (etype, nn, ecoord, gausses, section, stiffx, tt, t0, tdisp, rnqm )
1395 
1396  gausses(1)%nqm(1:12) = rnqm(1:12)
1397 
1398 ! write (6,'(a5,6a15)') 'dis-ij','x','y','z','theta-x','theta-y','theta-z'
1399 ! write (6,'(a,1p,6e15.5,0p)') 'dis-i',(tdisp(j),j= 1, 3),(tdisp(j),j= 7, 9)
1400 ! write (6,'(a,1p,6e15.5,0p)') 'dis-j',(tdisp(j),j= 4, 6),(tdisp(j),j=10,12)
1401 ! write (6,'(a5,6a15)') 'nqm-ij','N','Qy','QZ','Mx','My','Mz'
1402 ! write (6,'(a,1p,6e15.5,0p)') 'nqm-i',(rnqm(j),j= 1, 3),(rnqm(j),j= 7, 9)
1403 ! write (6,'(a,1p,6e15.5,0p)') 'nqm-j',(rnqm(j),j= 4, 6),(rnqm(j),j=10,12)
1404 ! write (6,'(a)') ''
1405 
1406  !--------------------------------------------------------------------
1407 
1408  return
1409 
1410  !####################################################################
1411  end subroutine nodalstress_beam_641
1412  !####################################################################
1413  ! > (Gaku Hashimoto, The University of Tokyo, 2013/09/13)
1414 
1415  !####################################################################
1416  subroutine elementalstress_beam_641 &
1417  ( gausses, estrain, estress, enqm )
1418  !####################################################################
1419  use m_fstr
1420  use mmechgauss
1421  implicit none
1422 
1423  !--------------------------------------------------------------------
1424 
1425  type(tgaussstatus), intent(inout) :: gausses(:)
1426  real(kind = kreal), intent(out) :: estrain(6)
1427  real(kind = kreal), intent(out) :: estress(6)
1428  real(kind = kreal), intent(out) :: enqm(12)
1429 
1430  !--------------------------------------------------------------------
1431 
1432  estrain(1:6) = gausses(1)%strain_out(1:6)
1433  estress(1:6) = gausses(1)%stress_out(1:6)
1434  enqm(1:12) = gausses(1)%nqm(1:12)
1435 
1436  end subroutine elementalstress_beam_641
1437 
1438  !####################################################################
1439  subroutine updatest_beam_641 &
1440  (etype, nn, ecoord, u, du, gausses, section, qf, tt, t0)
1441  !####################################################################
1442 
1443  use mmechgauss
1444 
1445  !--------------------------------------------------------------------
1446 
1447  integer, intent(in) :: etype
1448  integer, intent(in) :: nn
1449  real(kind=kreal), intent(in) :: ecoord(3, nn)
1450  real(kind=kreal), intent(in) :: u(3, nn)
1451  real(kind=kreal), intent(in) :: du(3, nn)
1452  type(tgaussstatus), intent(in) :: gausses(:)
1453  real(kind=kreal), intent(in) :: section(:)
1454  real(kind=kreal), intent(out) :: qf(nn*3)
1455  real(kind=kreal), intent(in), optional :: tt(nn), t0(nn)
1456 
1457  !--------------------------------------------------------------------
1458  real(kind = kreal) :: stiff(nn*3, nn*3), totaldisp(nn*3)
1459  integer(kind = kint) :: i
1460 
1461  call stf_beam_641(etype, nn, ecoord, gausses, section, stiff, tt, t0)
1462 
1463  totaldisp = 0.d0
1464  do i=1,nn
1465  totaldisp(3*i-2:3*i) = u(1:3,i) + du(1:3,i)
1466  end do
1467 
1468  qf = matmul(stiff,totaldisp)
1469 
1470  end subroutine updatest_beam_641
1471 
1472 end module
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
elementinfo::getnumberofnodes
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
Definition: element.f90:132
elementinfo::getsubface
subroutine getsubface(intype, innumber, outtype, nodes)
Find the definition of surface of the element.
Definition: element.f90:194
hecmw
Definition: hecmw.f90:6
m_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.F90:15
m_fstr::kel611timoshenko
integer(kind=kint), parameter kel611timoshenko
Definition: m_fstr.F90:87
m_fstr::ref_temp
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.F90:143
m_static_lib_beam
This module provide common functions of beam elements.
Definition: static_LIB_beam.f90:6
m_static_lib_beam::updatest_beam
subroutine updatest_beam(etype, nn, ecoord, u, du, section, gausses, QF, formulation)
Definition: static_LIB_beam.f90:184
m_static_lib_beam::nqm_beam_641
subroutine nqm_beam_641(etype, nn, ecoord, gausses, section, stiff, tt, t0, tdisp, rnqm)
Calculate N,Q,M vector of BEAM elements.
Definition: static_LIB_beam.f90:605
m_static_lib_beam::updatest_beam_641
subroutine updatest_beam_641(etype, nn, ecoord, u, du, gausses, section, qf, tt, t0)
Definition: static_LIB_beam.f90:1441
m_static_lib_beam::stf_beam_local
subroutine stf_beam_local(le, section, E, P, formulation, stiff)
Calculate the local stiffness matrix of 611 beam elements.
Definition: static_LIB_beam.f90:61
m_static_lib_beam::stf_beam
subroutine stf_beam(etype, nn, ecoord, section, E, P, STIFF, formulation)
Calculate stiff matrix of BEAM elements.
Definition: static_LIB_beam.f90:154
m_static_lib_beam::nqm_beam
subroutine nqm_beam(nn, ecoord, gausses, section, ul, rnqm, formulation)
Calculate elemental section force (N, Q, M) of 611 beam elements.
Definition: static_LIB_beam.f90:266
m_static_lib_beam::framtr
subroutine framtr(refx, xl, le, t)
Definition: static_LIB_beam.f90:20
m_static_lib_beam::elementalstress_beam
subroutine elementalstress_beam(gausses, estrain, estress, enqm)
Copy elemental stress/strain/NQM of 611 beam elements from Gauss status.
Definition: static_LIB_beam.f90:374
m_static_lib_beam::dl_beam_641
subroutine dl_beam_641(etype, nn, xx, yy, zz, rho, ltype, params, section, vect, nsize)
Definition: static_LIB_beam.f90:819
m_static_lib_beam::elementalstress_beam_641
subroutine elementalstress_beam_641(gausses, estrain, estress, enqm)
Definition: static_LIB_beam.f90:1418
m_static_lib_beam::updatest_beam_641_from_611
subroutine updatest_beam_641_from_611(ecoord, u, du, gausses, section, qf, formulation)
Definition: static_LIB_beam.f90:239
m_static_lib_beam::nodalstress_beam_641
subroutine nodalstress_beam_641(etype, nn, ecoord, gausses, section, edisp, ndstrain, ndstress, tt, t0, ntemp)
Definition: static_LIB_beam.f90:1083
m_static_lib_beam::tload_beam_641
subroutine tload_beam_641(etype, nn, ndof, xx, yy, zz, tt, t0, gausses, section, vect)
Definition: static_LIB_beam.f90:936
m_static_lib_beam::stf_beam_641_from_611
subroutine stf_beam_641_from_611(ecoord, gausses, section, stiff, formulation)
Definition: static_LIB_beam.f90:214
m_static_lib_beam::stf_beam_641
subroutine stf_beam_641(etype, nn, ecoord, gausses, section, stiff, tt, t0)
Calculate stiff matrix of BEAM elements.
Definition: static_LIB_beam.f90:391
m_static_lib_beam::nodalstress_beam
subroutine nodalstress_beam(etype, nn, ecoord, gausses, section, edisp, ndstrain, ndstress, formulation)
Calculate NODAL STRESS and STRAIN of 611 beam elements.
Definition: static_LIB_beam.f90:294
m_utilities
This module provides aux functions.
Definition: utilities.f90:6
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
mmechgauss::tgaussstatus
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:15