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