FrontISTR  6.0.0
Large-scale structural analysis program with finit element method
static_LIB_C3D8.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 !-------------------------------------------------------------------------------
7 module m_static_lib_c3d8
8 
9  use hecmw, only : kint, kreal
10  use elementinfo
11 
12  implicit none
13 
14 contains
15 
16 
22  !----------------------------------------------------------------------*
23  subroutine stf_c3d8bbar &
24  (etype, nn, ecoord, gausses, stiff, cdsys_id, coords, &
25  time, tincr, u, temperature)
26  !----------------------------------------------------------------------*
27 
28  use mmechgauss
29  use m_matmatrix
30  use m_common_struct
31  use m_static_lib_3d, only: geomat_c3
32 
33  !---------------------------------------------------------------------
34 
35  integer(kind=kint), intent(in) :: etype
36  integer(kind=kint), intent(in) :: nn
37  real(kind=kreal), intent(in) :: ecoord(3, nn)
38  type(tgaussstatus), intent(in) :: gausses(:)
39  real(kind=kreal), intent(out) :: stiff(:,:)
40  integer(kind=kint), intent(in) :: cdsys_id
41  real(kind=kreal), intent(inout) :: coords(3, 3)
42  real(kind=kreal), intent(in) :: time
43  real(kind=kreal), intent(in) :: tincr
44  real(kind=kreal), intent(in), optional :: u(:, :)
45  real(kind=kreal), intent(in) :: temperature(nn)
46 
47  !---------------------------------------------------------------------
48 
49  integer(kind=kint) :: flag
50  integer(kind=kint), parameter :: ndof = 3
51  real(kind=kreal) :: d(6, 6),b(6, ndof*nn),db(6, ndof*nn)
52  real(kind=kreal) :: gderiv(nn, 3),stress(6),mat(6, 6)
53  real(kind=kreal) :: det, wg, temp, spfunc(nn)
54  integer(kind=kint) :: i, j, lx, serr
55  real(kind=kreal) :: naturalcoord(3)
56  real(kind=kreal) :: gdispderiv(3, 3)
57  real(kind=kreal) :: b1(6, ndof*nn), bbar(nn, 3)
58  real(kind=kreal) :: smat(9, 9), elem(3, nn)
59  real(kind=kreal) :: bn(9, ndof*nn), sbn(9, ndof*nn)
60  real(kind=kreal) :: b4, b6, b8, coordsys(3, 3)
61 
62  !---------------------------------------------------------------------
63 
64  stiff(:, :) = 0.0d0
65  ! we suppose the same material type in the element
66  flag = gausses(1)%pMaterial%nlgeom_flag
67  if( .not. present(u) ) flag = infinitesimal ! enforce to infinitesimal deformation analysis
68  elem(:, :) = ecoord(:, :)
69  if( flag == updatelag ) elem(:, :) = ecoord(:, :)+u(:, :)
70 
71  ! dilatation component at centroid
72  naturalcoord = 0.0d0
73  call getglobalderiv(etype, nn, naturalcoord, elem, det, bbar)
74 
75  do lx = 1, numofquadpoints(etype)
76 
77  call getquadpoint( etype, lx, naturalcoord(:) )
78  call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
79 
80  if( cdsys_id > 0 ) then
81  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
82  if( serr == -1 ) stop "Fail to setup local coordinate"
83  if( serr == -2 ) then
84  !$omp critical
85  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
86  !$omp end critical
87  end if
88  end if
89 
90  call getshapefunc( etype, naturalcoord, spfunc )
91  temp = dot_product( temperature, spfunc )
92  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
93 
94  if( flag == updatelag ) then
95  call geomat_c3( gausses(lx)%stress, mat )
96  d(:, :) = d(:, :)-mat
97  endif
98 
99  wg = getweight(etype, lx)*det
100  b(1:6, 1:nn*ndof) = 0.0d0
101  do j = 1, nn
102  b4 = ( bbar(j, 1)-gderiv(j, 1) )/3.0d0
103  b6 = ( bbar(j, 2)-gderiv(j, 2) )/3.0d0
104  b8 = ( bbar(j, 3)-gderiv(j, 3) )/3.0d0
105  b(1, 3*j-2) = gderiv(j, 1)+b4
106  b(1, 3*j-1) = b6
107  b(1, 3*j ) = b8
108 
109  b(2, 3*j-2) = b4
110  b(2, 3*j-1) = gderiv(j, 2)+b6
111  b(2, 3*j ) = b8
112 
113  b(3, 3*j-2) = b4
114  b(3, 3*j-1) = b6
115  b(3, 3*j ) = gderiv(j, 3)+b8
116 
117  b(4, 3*j-2) = gderiv(j, 2)
118  b(4, 3*j-1) = gderiv(j, 1)
119  b(5, 3*j-1) = gderiv(j, 3)
120  b(5, 3*j ) = gderiv(j, 2)
121  b(6, 3*j-2) = gderiv(j, 3)
122  b(6, 3*j ) = gderiv(j, 1)
123  end do
124 
125  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
126  if( flag == totallag ) then
127  ! ---dudx(i,j) ==> gdispderiv(i,j)
128  gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
129  b1(1:6, 1:nn*ndof) = 0.0d0
130  do j = 1, nn
131  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
132  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
133  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
134  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
135  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
136  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
137  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
138  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
139  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
140  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
141  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
142  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
143  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
144  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
145  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
146  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
147  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
148  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
149  end do
150  ! ---BL = BL0 + BL1
151  do j = 1, nn*ndof
152  b(:, j) = b(:, j)+b1(:, j)
153  end do
154  end if
155 
156  db(1:6, 1:nn*ndof) = matmul( d, b(1:6, 1:nn*ndof) )
157  do j=1,nn*ndof
158  do i=1,nn*ndof
159  stiff(i, j) = stiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
160  end do
161  end do
162 
163  ! calculate the initial stress matrix
164  if( flag == totallag .or. flag == updatelag ) then
165  stress(1:6) = gausses(lx)%stress
166  bn(1:9, 1:nn*ndof) = 0.0d0
167  do j = 1, nn
168  bn(1, 3*j-2) = gderiv(j, 1)
169  bn(2, 3*j-1) = gderiv(j, 1)
170  bn(3, 3*j ) = gderiv(j, 1)
171  bn(4, 3*j-2) = gderiv(j, 2)
172  bn(5, 3*j-1) = gderiv(j, 2)
173  bn(6, 3*j ) = gderiv(j, 2)
174  bn(7, 3*j-2) = gderiv(j, 3)
175  bn(8, 3*j-1) = gderiv(j, 3)
176  bn(9, 3*j ) = gderiv(j, 3)
177  end do
178  smat(:, :) = 0.0d0
179  do j= 1, 3
180  smat(j , j ) = stress(1)
181  smat(j , j+3) = stress(4)
182  smat(j , j+6) = stress(6)
183  smat(j+3, j ) = stress(4)
184  smat(j+3, j+3) = stress(2)
185  smat(j+3, j+6) = stress(5)
186  smat(j+6, j ) = stress(6)
187  smat(j+6, j+3) = stress(5)
188  smat(j+6, j+6) = stress(3)
189  end do
190  sbn(1:9, 1:nn*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:nn*ndof) )
191  do j=1,nn*ndof
192  do i=1,nn*ndof
193  stiff(i, j) = stiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
194  end do
195  end do
196  end if
197 
198  end do ! gauss roop
199 
200  end subroutine stf_c3d8bbar
201 
202 
204  !----------------------------------------------------------------------*
205  subroutine update_c3d8bbar &
206  (etype, nn, ecoord, u, du, cdsys_id, coords, &
207  qf ,gausses, iter, time, tincr, tt,t0, tn )
208  !----------------------------------------------------------------------*
209 
210  use m_fstr
211  use mmaterial
212  use mmechgauss
213  use m_matmatrix
214  use m_elastoplastic
215  use mhyperelastic
216  use m_utilities
217  use m_static_lib_3d
218 
219  !---------------------------------------------------------------------
220 
221  integer(kind=kint), intent(in) :: etype
222  integer(kind=kint), intent(in) :: nn
223  real(kind=kreal), intent(in) :: ecoord(3, nn)
224  real(kind=kreal), intent(in) :: u(3, nn)
225  real(kind=kreal), intent(in) :: du(3, nn)
226  integer(kind=kint), intent(in) :: cdsys_id
227  real(kind=kreal), intent(inout) :: coords(3, 3)
228  real(kind=kreal), intent(out) :: qf(nn*3)
229  type(tgaussstatus), intent(inout) :: gausses(:)
230  integer, intent(in) :: iter
231  real(kind=kreal), intent(in) :: time
232  real(kind=kreal), intent(in) :: tincr
233  real(kind=kreal), intent(in) :: tt(nn)
234  real(kind=kreal), intent(in) :: t0(nn)
235  real(kind=kreal), intent(in) :: tn(nn)
236 
237  !---------------------------------------------------------------------
238 
239  integer(kind=kint) :: flag
240  integer(kind=kint), parameter :: ndof = 3
241  real(kind=kreal) :: b(6, ndof*nn), b1(6, ndof*nn)
242  real(kind=kreal) :: gderiv(nn, 3), gderiv1(nn,3), gdispderiv(3, 3), f(3,3), det, det1, wg
243  integer(kind=kint) :: j, lx, mtype, serr
244  real(kind=kreal) :: naturalcoord(3), rot(3, 3), spfunc(nn)
245  real(kind=kreal) :: totaldisp(3, nn), elem(3, nn), elem1(3, nn), coordsys(3, 3)
246  real(kind=kreal) :: dstrain(6)
247  real(kind=kreal) :: dvol, vol0, bbar(nn, 3), derivdum(1:ndof, 1:ndof), bbar2(nn, 3)
248  real(kind=kreal) :: b4, b6, b8, ttc, tt0, ttn, alpo(3), ina(1), epsth(6)
249  logical :: ierr, matlaniso
250 
251  !---------------------------------------------------------------------
252 
253  qf(:) = 0.0d0
254  ! we suppose the same material type in the element
255  flag = gausses(1)%pMaterial%nlgeom_flag
256  elem(:, :) = ecoord(:, :)
257  totaldisp(:, :) = u(:, :)+du(:, :)
258  if( flag == updatelag ) then
259  elem(:, :) = ( 0.5d0*du(:, :)+u(:, :) )+ecoord(:, :)
260  elem1(:, :) = ( du(:, :)+u(:, :) )+ecoord(:, :)
261  ! elem = elem1
262  totaldisp(:, :) = du(:, :)
263  end if
264 
265  matlaniso = .false.
266  ina = tt(1)
267  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
268  if( .not. ierr ) matlaniso = .true.
269 
270  ! dilatation at centroid
271  naturalcoord = 0.0d0
272  call getglobalderiv(etype, nn, naturalcoord, elem, det, bbar)
273  derivdum = matmul( totaldisp(1:ndof, 1:nn), bbar(1:nn, 1:ndof) )
274  vol0 = ( derivdum(1, 1)+derivdum(2, 2)+derivdum(3, 3) )/3.0d0
275  if( flag == updatelag ) call getglobalderiv(etype, nn, naturalcoord, elem1, det, bbar2)
276 
277  do lx = 1, numofquadpoints(etype)
278 
279  mtype = gausses(lx)%pMaterial%mtype
280 
281  call getquadpoint( etype, lx, naturalcoord(:) )
282  call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
283 
284  if( cdsys_id > 0 ) then
285  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
286  if( serr == -1 ) stop "Fail to setup local coordinate"
287  if( serr == -2 ) then
288  !$omp critical
289  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
290  !$omp end critical
291  end if
292  end if
293 
294  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
295  dvol = vol0-( gdispderiv(1, 1)+gdispderiv(2, 2)+gdispderiv(3, 3) )/3.0d0
296  !gdispderiv(1, 1) = gdispderiv(1, 1)+dvol
297  !gdispderiv(2, 2) = gdispderiv(2, 2)+dvol
298  !gdispderiv(3, 3) = gdispderiv(3, 3)+dvol
299 
300  ! ========================================================
301  ! UPDATE STRAIN and STRESS
302  ! ========================================================
303 
304  ! Thermal Strain
305  call getshapefunc(etype, naturalcoord, spfunc)
306  ttc = dot_product(tt, spfunc)
307  tt0 = dot_product(t0, spfunc)
308  ttn = dot_product(tn, spfunc)
309  call cal_thermal_expansion_c3( tt0, ttc, gausses(lx)%pMaterial, coordsys, matlaniso, epsth )
310 
311  ! Update strain
312  ! Small strain
313  dstrain(1) = gdispderiv(1, 1)+dvol
314  dstrain(2) = gdispderiv(2, 2)+dvol
315  dstrain(3) = gdispderiv(3, 3)+dvol
316  dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
317  dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
318  dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
319  dstrain(:) = dstrain(:)-epsth(:)
320 
321  f(1:3,1:3) = 0.d0; f(1,1)=1.d0; f(2,2)=1.d0; f(3,3)=1.d0; !deformation gradient
322  if( flag == infinitesimal ) then
323  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
324 
325  else if( flag == totallag ) then
326  ! Green-Lagrange strain
327  dstrain(1) = dstrain(1)+0.5d0*dot_product( gdispderiv(:, 1), gdispderiv(:, 1) )
328  dstrain(2) = dstrain(2)+0.5d0*dot_product( gdispderiv(:, 2), gdispderiv(:, 2) )
329  dstrain(3) = dstrain(3)+0.5d0*dot_product( gdispderiv(:, 3), gdispderiv(:, 3) )
330  dstrain(4) = dstrain(4)+( gdispderiv(1, 1)*gdispderiv(1, 2) &
331  +gdispderiv(2, 1)*gdispderiv(2, 2)+gdispderiv(3, 1)*gdispderiv(3, 2) )
332  dstrain(5) = dstrain(5)+( gdispderiv(1, 2)*gdispderiv(1, 3) &
333  +gdispderiv(2, 2)*gdispderiv(2, 3)+gdispderiv(3, 2)*gdispderiv(3, 3) )
334  dstrain(6) = dstrain(6)+( gdispderiv(1, 1)*gdispderiv(1, 3) &
335  +gdispderiv(2, 1)*gdispderiv(2, 3)+gdispderiv(3, 1)*gdispderiv(3, 3) )
336 
337  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
338  f(1:3,1:3) = f(1:3,1:3) + gdispderiv(1:3,1:3)
339 
340  else if( flag == updatelag ) then
341  rot = 0.0d0
342  rot(1, 2)= 0.5d0*( gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
343  rot(2, 3)= 0.5d0*( gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
344  rot(1, 3)= 0.5d0*( gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
345 
346  gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
347  call getglobalderiv(etype, nn, naturalcoord, ecoord, det1, gderiv1)
348  f(1:3,1:3) = f(1:3,1:3) + matmul( u(1:ndof, 1:nn)+du(1:ndof, 1:nn), gderiv1(1:nn, 1:ndof) )
349 
350  end if
351 
352  ! Update stress
353  call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn )
354 
355  ! ========================================================
356  ! calculate the internal force ( equivalent nodal force )
357  ! ========================================================
358  ! Small strain
359  b(1:6, 1:nn*ndof) = 0.0d0
360  do j = 1, nn
361  b4 = ( bbar(j, 1)-gderiv(j, 1) )/3.0d0
362  b6 = ( bbar(j, 2)-gderiv(j, 2) )/3.0d0
363  b8 = ( bbar(j, 3)-gderiv(j, 3) )/3.0d0
364  b(1, 3*j-2) = gderiv(j, 1)+b4
365  b(1, 3*j-1) = b6
366  b(1, 3*j ) = b8
367 
368  b(2, 3*j-2) = b4
369  b(2, 3*j-1) = gderiv(j, 2)+b6
370  b(2, 3*j ) = b8
371 
372  b(3, 3*j-2) = b4
373  b(3, 3*j-1) = b6
374  b(3, 3*j ) = gderiv(j, 3)+b8
375 
376  b(4, 3*j-2) = gderiv(j, 2)
377  b(4, 3*j-1) = gderiv(j, 1)
378  b(5, 3*j-1) = gderiv(j, 3)
379  b(5, 3*j ) = gderiv(j, 2)
380  b(6, 3*j-2) = gderiv(j, 3)
381  b(6, 3*j ) = gderiv(j, 1)
382  end do
383 
384  if( flag == infinitesimal ) then
385 
386  else if( flag == totallag ) then
387 
388  b1(1:6, 1:nn*ndof) = 0.0d0
389  do j = 1, nn
390  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
391  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
392  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
393  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
394  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
395  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
396  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
397  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
398  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
399  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
400  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
401  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
402  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
403  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
404  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
405  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
406  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
407  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
408  end do
409  ! BL = BL0 + BL1
410  do j = 1, nn*ndof
411  b(:, j) = b(:, j)+b1(:, j)
412  end do
413 
414  else if( flag == updatelag ) then
415 
416  call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv)
417  b(1:6, 1:nn*ndof) = 0.0d0
418  do j = 1, nn
419  b4 = ( bbar2(j, 1)-gderiv(j, 1) )/3.0d0
420  b6 = ( bbar2(j, 2)-gderiv(j, 2) )/3.0d0
421  b8 = ( bbar2(j, 3)-gderiv(j, 3) )/3.0d0
422  b(1, 3*j-2) = gderiv(j, 1)+b4
423  b(1, 3*j-1) = b6
424  b(1, 3*j ) = b8
425 
426  b(2, 3*j-2) = b4
427  b(2, 3*j-1) = gderiv(j, 2)+b6
428  b(2, 3*j ) = b8
429 
430  b(3, 3*j-2) = b4
431  b(3, 3*j-1) = b6
432  b(3, 3*j ) = gderiv(j, 3)+b8
433 
434  b(4, 3*j-2) = gderiv(j, 2)
435  b(4, 3*j-1) = gderiv(j, 1)
436  b(5, 3*j-1) = gderiv(j, 3)
437  b(5, 3*j ) = gderiv(j, 2)
438  b(6, 3*j-2) = gderiv(j, 3)
439  b(6, 3*j ) = gderiv(j, 1)
440  end do
441  end if
442 
443  !! calculate the Internal Force
444  wg = getweight(etype, lx)*det
445  qf(1:nn*ndof) = qf(1:nn*ndof) &
446  +matmul( gausses(lx)%stress(1:6), b(1:6, 1:nn*ndof) )*wg
447 
448  ! calculate strain energy
449  gausses(lx)%strain_energy = 0.5d0*dot_product(gausses(lx)%stress(1:6)+gausses(lx)%stress_bak(1:6), &
450  & gausses(lx)%strain(1:6)-gausses(lx)%strain_bak(1:6))*wg
451  end do
452 
453  end subroutine update_c3d8bbar
454 
456  !----------------------------------------------------------------------*
457  subroutine tload_c3d8bbar &
458  (etype, nn, xx, yy, zz, tt, t0, &
459  gausses, vect, cdsys_id, coords)
460  !----------------------------------------------------------------------*
461 
462  use m_fstr
463  use mmechgauss
464  use m_matmatrix
465  use m_utilities
466 
467  !---------------------------------------------------------------------
468 
469  integer(kind=kint), parameter :: ndof = 3
470  integer(kind=kint), intent(in) :: etype, nn
471  type(tgaussstatus), intent(in) :: gausses(:)
472  real(kind=kreal), intent(in) :: xx(nn), yy(nn), zz(nn)
473  real(kind=kreal), intent(in) :: tt(nn), t0(nn)
474  real(kind=kreal), intent(out) :: vect(nn*ndof)
475  integer(kind=kint), intent(in) :: cdsys_ID
476  real(kind=kreal), intent(inout) :: coords(3, 3)
477 
478  !---------------------------------------------------------------------
479 
480  real(kind=kreal) :: alp, alp0, d(6, 6), b(6, ndof*nn)
481  real(kind=kreal) :: b4, b6, b8, det, ecoord(3, nn)
482  integer(kind=kint) :: j, LX, serr
483  real(kind=kreal) :: estrain(6), sgm(6), h(nn)
484  real(kind=kreal) :: naturalcoord(3), gderiv(nn, 3)
485  real(kind=kreal) :: wg, outa(1), ina(1), bbar(nn, 3), alpo(3), alpo0(3), coordsys(3, 3)
486  real(kind=kreal) :: tempc, temp0, tempc0, temp00, thermal_eps, tm(6,6)
487  logical :: ierr, matlaniso
488 
489  !---------------------------------------------------------------------
490 
491  matlaniso = .false.
492 
493  if( cdsys_id > 0 ) then ! cannot define aniso expansion when no local coord defined
494  ina = tt(1)
495  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
496  if( .not. ierr ) matlaniso = .true.
497  end if
498 
499  vect(:) = 0.0d0
500 
501  ecoord(1, :) = xx(:)
502  ecoord(2, :) = yy(:)
503  ecoord(3, :) = zz(:)
504 
505  naturalcoord = 0.0d0
506  call getglobalderiv(etype, nn, naturalcoord, ecoord, det, bbar)
507  call getshapefunc( etype, naturalcoord, h(1:nn) )
508  tempc0 = dot_product( h(1:nn), tt(1:nn) )
509  temp00 = dot_product( h(1:nn), t0(1:nn) )
510 
511  ! LOOP FOR INTEGRATION POINTS
512  do lx = 1, numofquadpoints(etype)
513 
514  call getquadpoint( etype, lx, naturalcoord(:) )
515  call getshapefunc( etype, naturalcoord, h(1:nn) )
516  call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
517 
518  if( cdsys_id > 0 ) then
519  call set_localcoordsys(coords, g_localcoordsys(cdsys_id), coordsys, serr)
520  if( serr == -1 ) stop "Fail to setup local coordinate"
521  if( serr == -2 ) then
522  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
523  end if
524  end if
525 
526  ! WEIGHT VALUE AT GAUSSIAN POINT
527  wg = getweight(etype, lx)*det
528  b(1:6, 1:nn*ndof) = 0.0d0
529  do j = 1, nn
530  b4 = ( bbar(j, 1)-gderiv(j, 1) )/3.0d0
531  b6 = ( bbar(j, 2)-gderiv(j, 2) )/3.0d0
532  b8 = ( bbar(j, 3)-gderiv(j, 3) )/3.0d0
533  b(1, 3*j-2) = gderiv(j, 1)+b4
534  b(1, 3*j-1) = b6
535  b(1, 3*j ) = b8
536 
537  b(2, 3*j-2) = b4
538  b(2, 3*j-1) = gderiv(j, 2)+b6
539  b(2, 3*j ) = b8
540 
541  b(3, 3*j-2) = b4
542  b(3, 3*j-1) = b6
543  b(3, 3*j ) = gderiv(j, 3)+b8
544 
545  b(4, 3*j-2) = gderiv(j, 2)
546  b(4, 3*j-1) = gderiv(j, 1)
547  b(5, 3*j-1) = gderiv(j, 3)
548  b(5, 3*j ) = gderiv(j, 2)
549  b(6, 3*j-2) = gderiv(j, 3)
550  b(6, 3*j ) = gderiv(j, 1)
551  end do
552 
553  tempc = dot_product( h(1:nn),tt(1:nn) )
554  temp0 = dot_product( h(1:nn),t0(1:nn) )
555 
556  call matlmatrix( gausses(lx), d3, d, 1.d0, 0.0d0, coordsys, tempc )
557 
558  ina(1) = tempc
559  if( matlaniso ) then
560  call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo(:), ierr, ina )
561  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
562  else
563  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
564  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_expansion)
565  alp = outa(1)
566  end if
567  ina(1) = temp0
568  if( matlaniso ) then
569  call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo0(:), ierr, ina )
570  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
571  else
572  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
573  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_expansion)
574  alp0 = outa(1)
575  end if
576 
577  !**
578  !** THERMAL strain
579  !**
580  if( matlaniso ) then
581  do j = 1, 3
582  estrain(j) = alpo(j)*(tempc0-ref_temp)-alpo0(j)*(temp00-ref_temp)
583  end do
584  estrain(4:6) = 0.0d0
585  call transformation(coordsys, tm)
586  estrain(:) = matmul( estrain(:), tm ) ! to global coord
587  estrain(4:6) = estrain(4:6)*2.0d0
588  else
589  thermal_eps = alp*(tempc0-ref_temp)-alp0*(temp00-ref_temp)
590  estrain(1:3) = thermal_eps
591  estrain(4:6) = 0.0d0
592  end if
593 
594  !**
595  !** SET SGM {s}=[D]{e}
596  !**
597  sgm(:) = matmul( d(:, :), estrain(:) )
598 
599  !**
600  !** CALCULATE LOAD {F}=[B]T{e}
601  !**
602  vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(:), b(:, :) )*wg
603 
604  end do
605 
606  end subroutine tload_c3d8bbar
607 
608 
609 end module m_static_lib_c3d8
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
elementinfo::getshapefunc
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coordinate system.
Definition: element.f90:754
elementinfo::getquadpoint
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:539
elementinfo::getglobalderiv
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:848
elementinfo::getweight
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:585
elementinfo::numofquadpoints
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:451
hecmw
Definition: hecmw.f90:6
m_common_struct
This modules defines common structures for fem analysis.
Definition: m_common_struct.f90:6
m_common_struct::g_localcoordsys
type(tlocalcoordsys), dimension(:), pointer, save g_localcoordsys
Definition: m_common_struct.f90:19
m_common_struct::set_localcoordsys
subroutine set_localcoordsys(coords, coordsys, outsys, ierr)
setup of coordinate system
Definition: m_common_struct.f90:47
m_elastoplastic
This module provide functions for elastoplastic calculation.
Definition: Elastoplastic.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:147
m_matmatrix
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
m_matmatrix::matlmatrix
subroutine matlmatrix(gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp, hdflag)
Calculate constituive matrix.
Definition: calMatMatrix.f90:30
m_static_lib_3d
This module provides common functions of Solid elements.
Definition: static_LIB_3d.f90:6
m_static_lib_3d::geomat_c3
subroutine geomat_c3(stress, mat)
Definition: static_LIB_3d.f90:19
m_static_lib_3d::cal_thermal_expansion_c3
subroutine cal_thermal_expansion_c3(tt0, ttc, material, coordsys, matlaniso, EPSTH)
Definition: static_LIB_3d.f90:521
m_static_lib_3d::update_stress3d
subroutine update_stress3d(flag, gauss, rot, dstrain, F, coordsys, time, tincr, ttc, tt0, ttn, hdflag)
Definition: static_LIB_3d.f90:589
m_static_lib_c3d8
This module contains several strategy to free locking problem in Eight-node hexagonal element.
Definition: static_LIB_C3D8.f90:7
m_static_lib_c3d8::tload_c3d8bbar
subroutine tload_c3d8bbar(etype, nn, XX, YY, ZZ, TT, T0, gausses, VECT, cdsys_ID, coords)
This subroutien calculate thermal loading.
Definition: static_LIB_C3D8.f90:460
m_static_lib_c3d8::update_c3d8bbar
subroutine update_c3d8bbar(etype, nn, ecoord, u, du, cdsys_ID, coords, qf, gausses, iter, time, tincr, TT, T0, TN)
Update Strain stress of this element.
Definition: static_LIB_C3D8.f90:208
m_static_lib_c3d8::stf_c3d8bbar
subroutine stf_c3d8bbar(etype, nn, ecoord, gausses, stiff, cdsys_ID, coords, time, tincr, u, temperature)
This subroutine calculate stiff matrix using b-bar method.
Definition: static_LIB_C3D8.f90:26
m_utilities
This module provides aux functions.
Definition: utilities.f90:6
m_utilities::transformation
subroutine transformation(jacob, tm)
Definition: utilities.f90:392
mhyperelastic
This module provides functions for hyperelastic calculation.
Definition: Hyperelastic.f90:6
mmaterial
This module summarizes all information of material properties.
Definition: material.f90:6
mmaterial::totallag
integer(kind=kint), parameter totallag
Definition: material.f90:14
mmaterial::mc_orthoexp
character(len=dict_key_length) mc_orthoexp
Definition: material.f90:146
mmaterial::infinitesimal
integer(kind=kint), parameter infinitesimal
Definition: material.f90:13
mmaterial::updatelag
integer(kind=kint), parameter updatelag
Definition: material.f90:15
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