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