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