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