FrontISTR  5.7.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 
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  use gauss_integration
739 
740  !--------------------------------------------------------------------
741 
742  integer(kind = kint), intent(in) :: etype
743  integer(kind = kint), intent(in) :: nn
744  integer(kind = kint), intent(in) :: ndof
745  type(tgaussstatus), intent(in) :: gausses(:)
746  real(kind = kreal), intent(in) :: section(:)
747  real(kind = kreal), intent(in) :: xx(nn), yy(nn), zz(nn)
748  real(kind = kreal), intent(in) :: tt(nn), t0(nn)
749  real(kind = kreal), intent(out) :: vect(nn*ndof)
750 
751  !--------------------------------------------------------------------
752 
753  real(kind = kreal) :: tempc, temp0
754  real(kind = kreal) :: ecoord(3, nn)
755  real(kind = kreal) :: ec(3, 2)
756  real(kind = kreal) :: ina1(1), outa1(2)
757  real(kind = kreal) :: ina2(1), outa2(1)
758  real(kind = kreal) :: alp, alp0
759  real(kind = kreal) :: ee, pp
760  real(kind = kreal) :: a
761  real(kind = kreal) :: refv(3)
762  real(kind = kreal) :: g
763  real(kind = kreal) :: le
764  real(kind = kreal) :: trans(3, 3), transt(3, 3)
765 
766  logical :: ierr
767 
768  !--------------------------------------------------------------------
769 
770  ecoord(1, 1:nn) = xx(1:nn)
771  ecoord(2, 1:nn) = yy(1:nn)
772  ecoord(3, 1:nn) = zz(1:nn)
773 
774  !--------------------------------------------------------------------
775 
776  tempc = 0.5d0*( tt(1)+tt(2) )
777  temp0 = 0.5d0*( t0(1)+t0(2) )
778 
779  !--------------------------------------------------------------
780 
781  ina1(1) = tempc
782 
783  call fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa1, ierr, ina1 )
784 
785  if( ierr ) then
786 
787  ee = gausses(1)%pMaterial%variables(m_youngs)
788  pp = gausses(1)%pMaterial%variables(m_poisson)
789 
790  else
791 
792  ee = outa1(1)
793  pp = outa1(2)
794 
795  end if
796 
797  !--------------------------------------------------------------
798 
799  ina2(1) = tempc
800 
801  call fetch_tabledata( mc_themoexp, gausses(1)%pMaterial%dict, outa2(:), ierr, ina2 )
802 
803  if( ierr ) stop "Fails in fetching expansion coefficient!"
804 
805  alp = outa2(1)
806 
807  !--------------------------------------------------------------
808 
809  ina2(1) = temp0
810 
811  call fetch_tabledata( mc_themoexp, gausses(1)%pMaterial%dict, outa2(:), ierr, ina2 )
812 
813  if( ierr ) stop "Fails in fetching expansion coefficient!"
814 
815  alp0 = outa2(1)
816 
817  !--------------------------------------------------------------------
818 
819  refv(1) = section(1)
820  refv(2) = section(2)
821  refv(3) = section(3)
822 
823  ec(1, 1) = ecoord(1, 1)
824  ec(2, 1) = ecoord(2, 1)
825  ec(3, 1) = ecoord(3, 1)
826  ec(1, 2) = ecoord(1, 2)
827  ec(2, 2) = ecoord(2, 2)
828  ec(3, 2) = ecoord(3, 2)
829 
830  call framtr(refv, ec, le, trans)
831 
832  transt= transpose( trans )
833 
834  !--------------------------------------------------------------------
835 
836  a = section(4)
837 
838  g = ee/( 2.0d0*( 1.0d0+pp ))
839 
840  !--------------------------------------------------------------------
841 
842  vect( 1) = -a*ee*( alp*( tempc-ref_temp )-alp0*( temp0-ref_temp ) )
843  vect( 2) = 0.0d0
844  vect( 3) = 0.0d0
845 
846  vect( 4) = a*ee*( alp*( tempc-ref_temp )-alp0*( temp0-ref_temp ) )
847  vect( 5) = 0.0d0
848  vect( 6) = 0.0d0
849 
850  vect( 7) = 0.0d0
851  vect( 8) = 0.0d0
852  vect( 9) = 0.0d0
853 
854  vect(10) = 0.0d0
855  vect(11) = 0.0d0
856  vect(12) = 0.0d0
857 
858  !--------------------------------------------------------------------
859 
860  vect( 1:3) = matmul( transt, vect(1:3) )
861  vect( 4:6) = matmul( transt, vect(4:6) )
862  vect( 7:9) = matmul( transt, vect(7:9) )
863  vect(10:12) = matmul( transt, vect(10:12) )
864 
865  !--------------------------------------------------------------------
866 
867  return
868 
869  !####################################################################
870  end subroutine tload_beam_641
871  !####################################################################
872  ! > (Gaku Hashimoto, The University of Tokyo, 2013/09/13)
873 
874 
875  ! (Gaku Hashimoto, The University of Tokyo, 2014/02/06) <
876  !####################################################################
877  subroutine nodalstress_beam_641 &
878  (etype, nn, ecoord, gausses, section, edisp, &
879  ndstrain, ndstress, tt, t0, ntemp)
880  !####################################################################
881 
882  use m_fstr
883  use mmechgauss
884 
885  !--------------------------------------------------------------------
886 
887  integer(kind = kint), intent(in) :: etype
888  integer(kind = kint), intent(in) :: nn
889  real(kind = kreal), intent(in) :: ecoord(3, nn)
890  type(tgaussstatus), intent(inout) :: gausses(:)
891  real(kind = kreal), intent(in) :: section(:)
892  real(kind = kreal), intent(in) :: edisp(3, nn)
893  real(kind = kreal), intent(out) :: ndstrain(nn, 6)
894  real(kind = kreal), intent(out) :: ndstress(nn, 6)
895  real(kind=kreal), intent(in), optional :: tt(nn), t0(nn)
896  integer(kind = kint), intent(in) :: ntemp
897 
898  !--------------------------------------------------------------------
899 
900  real(kind=kreal) :: stiffx(12, 12)
901  real(kind=kreal) :: tdisp(12)
902  real(kind=kreal) :: rnqm(12)
903 
904  !--------------------------------------------------------------------
905 
906  integer(kind = kint) :: i, j, k, jj
907 
908  real(kind = kreal) :: tempc, temp0
909  real(kind = kreal) :: ina1(1), outa1(2)
910  real(kind = kreal) :: ina2(1), outa2(1)
911  real(kind = kreal) :: alp, alp0
912  real(kind = kreal) :: ee, pp
913  real(kind = kreal) :: a, radius, angle(6)
914  real(kind = kreal) :: refv(3)
915  real(kind = kreal) :: le, l2, l3
916  real(kind = kreal) :: trans(3, 3), transt(3, 3)
917  real(kind = kreal) :: edisp_hat(3, nn)
918  real(kind = kreal) :: ec(3, 2)
919  real(kind = kreal) :: t(3, 3), t_hat(3, 3)
920  real(kind = kreal) :: t_hat_tmp(3, 3)
921  real(kind = kreal) :: e(3, 3), e_hat(3, 3)
922  real(kind = kreal) :: e_hat_tmp(3, 3)
923  real(kind = kreal) :: x1_hat, x2_hat, x3_hat
924  real(kind = kreal) :: pi
925 
926  logical :: ierr
927 
928  alp = 0.0d0; alp0 = 0.0d0
929  tempc = 0.0d0; temp0 = 0.0d0
930 
931  !--------------------------------------------------------------------
932 
933  pi = 4.0d0*datan( 1.0d0 )
934 
935  !--------------------------------------------------------------------
936 
937  if( present( tt ) .AND. present( t0 ) ) then
938 
939  tempc = 0.5d0*( tt(1)+tt(2) )
940  temp0 = 0.5d0*( t0(1)+t0(2) )
941 
942  end if
943 
944  !--------------------------------------------------------------------
945 
946  if( ntemp .EQ. 1 ) then
947 
948  ina1(1) = tempc
949 
950  call fetch_tabledata( mc_isoelastic, gausses(1)%pMaterial%dict, outa1, ierr, ina1 )
951 
952  else
953 
954  ierr = .true.
955 
956  end if
957 
958  !--------------------------------------------------------------
959 
960  if( ierr ) then
961 
962  ee = gausses(1)%pMaterial%variables(m_youngs)
963  pp = gausses(1)%pMaterial%variables(m_poisson)
964 
965  else
966 
967  ee = outa1(1)
968  pp = outa1(2)
969 
970  end if
971 
972  !--------------------------------------------------------------------
973 
974  if( ntemp .EQ. 1 ) then
975 
976  ina2(1) = tempc
977 
978  call fetch_tabledata( mc_themoexp, gausses(1)%pMaterial%dict, outa2(:), ierr, ina2 )
979 
980  if( ierr ) stop "Fails in fetching expansion coefficient!"
981 
982  alp = outa2(1)
983 
984  end if
985 
986  !--------------------------------------------------------------
987 
988  if( ntemp .EQ. 1 ) then
989 
990  ina2(1) = temp0
991 
992  call fetch_tabledata( mc_themoexp, gausses(1)%pMaterial%dict, outa2(:), ierr, ina2 )
993 
994  if( ierr ) stop "Fails in fetching expansion coefficient!"
995 
996  alp0 = outa2(1)
997 
998  end if
999 
1000  !--------------------------------------------------------------------
1001 
1002  refv(1) = section(1)
1003  refv(2) = section(2)
1004  refv(3) = section(3)
1005 
1006  ec(1, 1) = ecoord(1, 1)
1007  ec(2, 1) = ecoord(2, 1)
1008  ec(3, 1) = ecoord(3, 1)
1009  ec(1, 2) = ecoord(1, 2)
1010  ec(2, 2) = ecoord(2, 2)
1011  ec(3, 2) = ecoord(3, 2)
1012 
1013  call framtr(refv, ec, le, trans)
1014 
1015  transt= transpose( trans )
1016 
1017  l2 = le*le
1018  l3 = l2*le
1019 
1020  !--------------------------------------------------------------------
1021 
1022  a = section(4)
1023 
1024  radius = gausses(1)%pMaterial%variables(m_beam_radius)
1025 
1026  angle(1) = gausses(1)%pMaterial%variables(m_beam_angle1)
1027  angle(2) = gausses(1)%pMaterial%variables(m_beam_angle2)
1028  angle(3) = gausses(1)%pMaterial%variables(m_beam_angle3)
1029  angle(4) = gausses(1)%pMaterial%variables(m_beam_angle4)
1030  angle(5) = gausses(1)%pMaterial%variables(m_beam_angle5)
1031  angle(6) = gausses(1)%pMaterial%variables(m_beam_angle6)
1032 
1033  !--------------------------------------------------------------------
1034 
1035  do k = 1, 6
1036 
1037  !--------------------------------------------------------
1038 
1039  angle(k) = angle(k)/180.0d0*pi
1040 
1041  x2_hat = radius*dcos( angle(k) )
1042  x3_hat = radius*dsin( angle(k) )
1043 
1044  !--------------------------------------------------------
1045 
1046  jj = 0
1047  do j = 1, nn
1048 
1049  do i = 1, 3
1050 
1051  edisp_hat(i, j) = trans(i, 1)*edisp(1, j) &
1052  +trans(i, 2)*edisp(2, j) &
1053  +trans(i, 3)*edisp(3, j)
1054 
1055  jj = jj + 1
1056  tdisp(jj) = edisp(i,j)
1057 
1058  end do
1059 
1060  end do
1061 
1062  !--------------------------------------------------------
1063 
1064  x1_hat = 0.5d0*le
1065 
1066  e_hat = 0.0d0
1067  e_hat(1, 1) = ( edisp_hat(1, 2)-edisp_hat(1, 1) )/le
1068 
1069  t_hat = 0.0d0
1070  t_hat(1, 1) = ee*( edisp_hat(1, 2)-edisp_hat(1, 1) )/le &
1071  -ee*x2_hat*( ( -6.0d0/l2+12.0d0*x1_hat/l3 )*edisp_hat(2, 1) &
1072  +( -4.0d0/le+6.0d0*x1_hat/l2 )*edisp_hat(3, 3) &
1073  +( 6.0d0/l2-12.0d0*x1_hat/l3 )*edisp_hat(2, 2) &
1074  +( -2.0d0/le+6.0d0*x1_hat/l2 )*edisp_hat(3, 4) ) &
1075  -ee*x3_hat*( ( -6.0d0/l2+12.0d0*x1_hat/l3 )*edisp_hat(3, 1) &
1076  +( 4.0d0/le-6.0d0*x1_hat/l2 )*edisp_hat(2, 3) &
1077  +( 6.0d0/l2-12.0d0*x1_hat/l3 )*edisp_hat(3, 2) &
1078  +( 2.0d0/le-6.0d0*x1_hat/l2 )*edisp_hat(2, 4) )
1079 
1080  if( ntemp .EQ. 1 ) then
1081 
1082  t_hat(1, 1) &
1083  = t_hat(1, 1) &
1084  -ee*( alp*( tempc-ref_temp )-alp0*( temp0-ref_temp ) )
1085 
1086  end if
1087 
1088  e_hat_tmp(1:3,:) = matmul( trans, e_hat(1:3,:) )
1089  t_hat_tmp(1:3,:) = matmul( trans, t_hat(1:3,:) )
1090 
1091  e(:, 1:3) = matmul( e_hat_tmp(:,1:3), transt )
1092  t(:, 1:3) = matmul( t_hat_tmp(:,1:3), transt )
1093 
1094  gausses(1)%strain(k) = e_hat(1, 1)
1095  gausses(1)%stress(k) = t_hat(1, 1)
1096 
1097  !set stress and strain for output
1098  gausses(1)%strain_out(k) = gausses(1)%strain(k)
1099  gausses(1)%stress_out(k) = gausses(1)%stress(k)
1100 
1101  !--------------------------------------------------------
1102 
1103  ndstrain(1, k) = 0.0d0
1104  ndstrain(2, k) = 0.0d0
1105  ndstrain(3, k) = 0.0d0
1106  ndstrain(4, k) = 0.0d0
1107 
1108  ndstress(1, k) = 0.0d0
1109  ndstress(2, k) = 0.0d0
1110  ndstress(3, k) = 0.0d0
1111  ndstress(4, k) = 0.0d0
1112 
1113  !--------------------------------------------------------
1114 
1115  x1_hat = 0.0d0
1116 
1117  e_hat = 0.0d0
1118  e_hat(1, 1) = ( edisp_hat(1, 2)-edisp_hat(1, 1) )/le
1119 
1120  t_hat = 0.0d0
1121  t_hat(1, 1) = ee*( edisp_hat(1, 2)-edisp_hat(1, 1) )/le &
1122  -ee*x2_hat*( ( -6.0d0/l2+12.0d0*x1_hat/l3 )*edisp_hat(2, 1) &
1123  +( -4.0d0/le+6.0d0*x1_hat/l2 )*edisp_hat(3, 3) &
1124  +( 6.0d0/l2-12.0d0*x1_hat/l3 )*edisp_hat(2, 2) &
1125  +( -2.0d0/le+6.0d0*x1_hat/l2 )*edisp_hat(3, 4) ) &
1126  -ee*x3_hat*( ( -6.0d0/l2+12.0d0*x1_hat/l3 )*edisp_hat(3, 1) &
1127  +( 4.0d0/le-6.0d0*x1_hat/l2 )*edisp_hat(2, 3) &
1128  +( 6.0d0/l2-12.0d0*x1_hat/l3 )*edisp_hat(3, 2) &
1129  +( 2.0d0/le-6.0d0*x1_hat/l2 )*edisp_hat(2, 4) )
1130 
1131  if( ntemp .EQ. 1 ) then
1132 
1133  t_hat(1, 1) &
1134  = t_hat(1, 1) &
1135  -ee*( alp*( tempc-ref_temp )-alp0*( temp0-ref_temp ) )
1136 
1137  end if
1138 
1139  e_hat_tmp(1:3, :) = matmul( trans, e_hat(1:3, :) )
1140  t_hat_tmp(1:3, :) = matmul( trans, t_hat(1:3, :) )
1141 
1142  e(:, 1:3) = matmul( e_hat_tmp(:, 1:3), transt )
1143  t(:, 1:3) = matmul( t_hat_tmp(:, 1:3), transt )
1144 
1145  ndstrain(1, k) = e_hat(1, 1)
1146  ndstress(1, k) = t_hat(1, 1)
1147 
1148  !--------------------------------------------------------
1149 
1150  x1_hat = le
1151 
1152  e_hat = 0.0d0
1153  e_hat(1, 1) = ( edisp_hat(1, 2)-edisp_hat(1, 1) )/le
1154 
1155  t_hat = 0.0d0
1156  t_hat(1, 1) = ee*( edisp_hat(1, 2)-edisp_hat(1, 1) )/le &
1157  -ee*x2_hat*( ( -6.0d0/l2+12.0d0*x1_hat/l3 )*edisp_hat(2, 1) &
1158  +( -4.0d0/le+6.0d0*x1_hat/l2 )*edisp_hat(3, 3) &
1159  +( 6.0d0/l2-12.0d0*x1_hat/l3 )*edisp_hat(2, 2) &
1160  +( -2.0d0/le+6.0d0*x1_hat/l2 )*edisp_hat(3, 4) ) &
1161  -ee*x3_hat*( ( -6.0d0/l2+12.0d0*x1_hat/l3 )*edisp_hat(3, 1) &
1162  +( 4.0d0/le-6.0d0*x1_hat/l2 )*edisp_hat(2, 3) &
1163  +( 6.0d0/l2-12.0d0*x1_hat/l3 )*edisp_hat(3, 2) &
1164  +( 2.0d0/le-6.0d0*x1_hat/l2 )*edisp_hat(2, 4) )
1165 
1166  if( ntemp .EQ. 1 ) then
1167 
1168  t_hat(1, 1) &
1169  = t_hat(1, 1) &
1170  -ee*( alp*( tempc-ref_temp )-alp0*( temp0-ref_temp ) )
1171 
1172  end if
1173 
1174  e_hat_tmp(1:3, :) = matmul( trans, e_hat(1:3, :) )
1175  t_hat_tmp(1:3, :) = matmul( trans, t_hat(1:3, :) )
1176 
1177  e(:, 1:3) = matmul( e_hat_tmp(:, 1:3), transt )
1178  t(:, 1:3) = matmul( t_hat_tmp(:, 1:3), transt )
1179 
1180  ndstrain(2, k) = e_hat(1, 1)
1181  ndstress(2, k) = t_hat(1, 1)
1182 
1183  !--------------------------------------------------------
1184 
1185  end do
1186 
1187  !--------------------------------------------------------------------
1188  stiffx = 0.0
1189 
1190  call nqm_beam_641 &
1191  (etype, nn, ecoord, gausses, section, stiffx, tt, t0, tdisp, rnqm )
1192 
1193  gausses(1)%nqm(1:12) = rnqm(1:12)
1194 
1195 ! write (6,'(a5,6a15)') 'dis-ij','x','y','z','theta-x','theta-y','theta-z'
1196 ! write (6,'(a,1p,6e15.5,0p)') 'dis-i',(tdisp(j),j= 1, 3),(tdisp(j),j= 7, 9)
1197 ! write (6,'(a,1p,6e15.5,0p)') 'dis-j',(tdisp(j),j= 4, 6),(tdisp(j),j=10,12)
1198 ! write (6,'(a5,6a15)') 'nqm-ij','N','Qy','QZ','Mx','My','Mz'
1199 ! write (6,'(a,1p,6e15.5,0p)') 'nqm-i',(rnqm(j),j= 1, 3),(rnqm(j),j= 7, 9)
1200 ! write (6,'(a,1p,6e15.5,0p)') 'nqm-j',(rnqm(j),j= 4, 6),(rnqm(j),j=10,12)
1201 ! write (6,'(a)') ''
1202 
1203  !--------------------------------------------------------------------
1204 
1205  return
1206 
1207  !####################################################################
1208  end subroutine nodalstress_beam_641
1209  !####################################################################
1210  ! > (Gaku Hashimoto, The University of Tokyo, 2013/09/13)
1211 
1212  !####################################################################
1213  subroutine elementalstress_beam_641 &
1214  ( gausses, estrain, estress, enqm )
1215  !####################################################################
1216  use m_fstr
1217  use mmechgauss
1218  implicit none
1219 
1220  !--------------------------------------------------------------------
1221 
1222  type(tgaussstatus), intent(inout) :: gausses(:)
1223  real(kind = kreal), intent(out) :: estrain(6)
1224  real(kind = kreal), intent(out) :: estress(6)
1225  real(kind = kreal), intent(out) :: enqm(12)
1226 
1227  !--------------------------------------------------------------------
1228 
1229  estrain(1:6) = gausses(1)%strain_out(1:6)
1230  estress(1:6) = gausses(1)%stress_out(1:6)
1231  enqm(1:12) = gausses(1)%nqm(1:12)
1232 
1233  end subroutine elementalstress_beam_641
1234 
1235  !####################################################################
1236  subroutine updatest_beam_641 &
1237  (etype, nn, ecoord, u, du, gausses, section, qf, tt, t0)
1238  !####################################################################
1239 
1240  use mmechgauss
1241 
1242  !--------------------------------------------------------------------
1243 
1244  integer, intent(in) :: etype
1245  integer, intent(in) :: nn
1246  real(kind=kreal), intent(in) :: ecoord(3, nn)
1247  real(kind=kreal), intent(in) :: u(3, nn)
1248  real(kind=kreal), intent(in) :: du(3, nn)
1249  type(tgaussstatus), intent(in) :: gausses(:)
1250  real(kind=kreal), intent(in) :: section(:)
1251  real(kind=kreal), intent(out) :: qf(nn*3)
1252  real(kind=kreal), intent(in), optional :: tt(nn), t0(nn)
1253 
1254  !--------------------------------------------------------------------
1255  real(kind = kreal) :: stiff(nn*3, nn*3), totaldisp(nn*3)
1256  integer(kind = kint) :: i
1257 
1258  call stf_beam_641(etype, nn, ecoord, gausses, section, stiff, tt, t0)
1259 
1260  totaldisp = 0.d0
1261  do i=1,nn
1262  totaldisp(3*i-2:3*i) = u(1:3,i) + du(1:3,i)
1263  end do
1264 
1265  qf = matmul(stiff,totaldisp)
1266 
1267  end subroutine updatest_beam_641
1268 
1269 end module
m_utilities
This module provides aux functions.
Definition: utilities.f90:6
m_static_lib_beam
This module provide common functions of beam elements.
Definition: static_LIB_beam.f90:6
m_static_lib_beam::framtr
subroutine framtr(refx, xl, le, t)
Definition: static_LIB_beam.f90:19
gauss_integration
This module provides data for gauss quadrature.
Definition: GaussM.f90:6
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:880
m_static_lib_beam::elementalstress_beam_641
subroutine elementalstress_beam_641(gausses, estrain, estress, enqm)
Definition: static_LIB_beam.f90:1215
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
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_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::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
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
hecmw
Definition: hecmw.f90:6
mmechgauss::tgaussstatus
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13
elementinfo::getsubface
subroutine getsubface(intype, innumber, outtype, nodes)
Find the definition of surface of the element.
Definition: element.f90:193
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::updatest_beam
subroutine updatest_beam(etype, nn, ecoord, u, du, section, gausses, QF)
Definition: static_LIB_beam.f90:154
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:1238
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::dl_beam_641
subroutine dl_beam_641(etype, nn, xx, yy, zz, rho, ltype, params, section, vect, nsize)
Definition: static_LIB_beam.f90:615
elementinfo::getnumberofnodes
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
Definition: element.f90:131