FrontISTR  5.2.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 !-------------------------------------------------------------------------------
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), optional :: 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  if( present(temperature) ) then
89  call getshapefunc( etype, naturalcoord, spfunc )
90  temp = dot_product( temperature, spfunc )
91  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
92  else
93  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys )
94  end if
95 
96  if( flag == updatelag ) then
97  call geomat_c3( gausses(lx)%stress, mat )
98  d(:, :) = d(:, :)-mat
99  endif
100 
101  wg = getweight(etype, lx)*det
102  b(1:6, 1:nn*ndof) = 0.0d0
103  do j = 1, nn
104  b4 = ( bbar(j, 1)-gderiv(j, 1) )/3.0d0
105  b6 = ( bbar(j, 2)-gderiv(j, 2) )/3.0d0
106  b8 = ( bbar(j, 3)-gderiv(j, 3) )/3.0d0
107  b(1, 3*j-2) = gderiv(j, 1)+b4
108  b(1, 3*j-1) = b6
109  b(1, 3*j ) = b8
110 
111  b(2, 3*j-2) = b4
112  b(2, 3*j-1) = gderiv(j, 2)+b6
113  b(2, 3*j ) = b8
114 
115  b(3, 3*j-2) = b4
116  b(3, 3*j-1) = b6
117  b(3, 3*j ) = gderiv(j, 3)+b8
118 
119  b(4, 3*j-2) = gderiv(j, 2)
120  b(4, 3*j-1) = gderiv(j, 1)
121  b(5, 3*j-1) = gderiv(j, 3)
122  b(5, 3*j ) = gderiv(j, 2)
123  b(6, 3*j-2) = gderiv(j, 3)
124  b(6, 3*j ) = gderiv(j, 1)
125  end do
126 
127  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
128  if( flag == totallag ) then
129  ! ---dudx(i,j) ==> gdispderiv(i,j)
130  gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
131  b1(1:6, 1:nn*ndof) = 0.0d0
132  do j = 1, nn
133  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
134  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
135  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
136  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
137  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
138  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
139  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
140  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
141  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
142  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
143  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
144  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
145  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
146  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
147  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
148  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
149  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
150  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
151  end do
152  ! ---BL = BL0 + BL1
153  do j = 1, nn*ndof
154  b(:, j) = b(:, j)+b1(:, j)
155  end do
156  end if
157 
158  db(1:6, 1:nn*ndof) = matmul( d, b(1:6, 1:nn*ndof) )
159  forall( i=1:nn*ndof, j=1:nn*ndof )
160  stiff(i, j) = stiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
161  end forall
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  forall( i=1:nn*ndof, j=1:nn*ndof )
192  stiff(i, j) = stiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
193  end forall
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), optional :: tt(nn)
232  real(kind=kreal), intent(in), optional :: t0(nn)
233  real(kind=kreal), intent(in), optional :: 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  if( cdsys_id > 0 .AND. present(tt) ) then
265  ina = tt(1)
266  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
267  if( .not. ierr ) matlaniso = .true.
268  end if
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  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
289  end if
290  end if
291 
292  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
293  dvol = vol0-( gdispderiv(1, 1)+gdispderiv(2, 2)+gdispderiv(3, 3) )/3.0d0
294  !gdispderiv(1, 1) = gdispderiv(1, 1)+dvol
295  !gdispderiv(2, 2) = gdispderiv(2, 2)+dvol
296  !gdispderiv(3, 3) = gdispderiv(3, 3)+dvol
297 
298  ! ========================================================
299  ! UPDATE STRAIN and STRESS
300  ! ========================================================
301 
302  ! Thermal Strain
303  epsth = 0.0d0
304  if( present(tt) .AND. present(t0) ) then
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  end if
311 
312  ! Update strain
313  ! Small strain
314  dstrain(1) = gdispderiv(1, 1)+dvol
315  dstrain(2) = gdispderiv(2, 2)+dvol
316  dstrain(3) = gdispderiv(3, 3)+dvol
317  dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
318  dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
319  dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
320  dstrain(:) = dstrain(:)-epsth(:)
321 
322  f(1:3,1:3) = 0.d0; f(1,1)=1.d0; f(2,2)=1.d0; f(3,3)=1.d0; !deformation gradient
323  if( flag == infinitesimal ) then
324  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
325 
326  else if( flag == totallag ) then
327  ! Green-Lagrange strain
328  dstrain(1) = dstrain(1)+0.5d0*dot_product( gdispderiv(:, 1), gdispderiv(:, 1) )
329  dstrain(2) = dstrain(2)+0.5d0*dot_product( gdispderiv(:, 2), gdispderiv(:, 2) )
330  dstrain(3) = dstrain(3)+0.5d0*dot_product( gdispderiv(:, 3), gdispderiv(:, 3) )
331  dstrain(4) = dstrain(4)+( gdispderiv(1, 1)*gdispderiv(1, 2) &
332  +gdispderiv(2, 1)*gdispderiv(2, 2)+gdispderiv(3, 1)*gdispderiv(3, 2) )
333  dstrain(5) = dstrain(5)+( gdispderiv(1, 2)*gdispderiv(1, 3) &
334  +gdispderiv(2, 2)*gdispderiv(2, 3)+gdispderiv(3, 2)*gdispderiv(3, 3) )
335  dstrain(6) = dstrain(6)+( gdispderiv(1, 1)*gdispderiv(1, 3) &
336  +gdispderiv(2, 1)*gdispderiv(2, 3)+gdispderiv(3, 1)*gdispderiv(3, 3) )
337 
338  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
339  f(1:3,1:3) = f(1:3,1:3) + gdispderiv(1:3,1:3)
340 
341  else if( flag == updatelag ) then
342  rot = 0.0d0
343  rot(1, 2)= 0.5d0*( gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
344  rot(2, 3)= 0.5d0*( gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
345  rot(1, 3)= 0.5d0*( gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
346 
347  gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
348  call getglobalderiv(etype, nn, naturalcoord, ecoord, det1, gderiv1)
349  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) )
350 
351  end if
352 
353  ! Update stress
354  if( present(tt) .AND. present(t0) ) then
355  call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn )
356  else
357  call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr )
358  end if
359 
360  ! ========================================================
361  ! calculate the internal force ( equivalent nodal force )
362  ! ========================================================
363  ! Small strain
364  b(1:6, 1:nn*ndof) = 0.0d0
365  do j = 1, nn
366  b4 = ( bbar(j, 1)-gderiv(j, 1) )/3.0d0
367  b6 = ( bbar(j, 2)-gderiv(j, 2) )/3.0d0
368  b8 = ( bbar(j, 3)-gderiv(j, 3) )/3.0d0
369  b(1, 3*j-2) = gderiv(j, 1)+b4
370  b(1, 3*j-1) = b6
371  b(1, 3*j ) = b8
372 
373  b(2, 3*j-2) = b4
374  b(2, 3*j-1) = gderiv(j, 2)+b6
375  b(2, 3*j ) = b8
376 
377  b(3, 3*j-2) = b4
378  b(3, 3*j-1) = b6
379  b(3, 3*j ) = gderiv(j, 3)+b8
380 
381  b(4, 3*j-2) = gderiv(j, 2)
382  b(4, 3*j-1) = gderiv(j, 1)
383  b(5, 3*j-1) = gderiv(j, 3)
384  b(5, 3*j ) = gderiv(j, 2)
385  b(6, 3*j-2) = gderiv(j, 3)
386  b(6, 3*j ) = gderiv(j, 1)
387  end do
388 
389  if( flag == infinitesimal ) then
390 
391  else if( flag == totallag ) then
392 
393  b1(1:6, 1:nn*ndof) = 0.0d0
394  do j = 1, nn
395  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
396  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
397  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
398  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
399  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
400  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
401  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
402  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
403  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
404  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
405  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
406  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
407  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
408  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
409  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
410  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
411  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
412  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
413  end do
414  ! BL = BL0 + BL1
415  do j = 1, nn*ndof
416  b(:, j) = b(:, j)+b1(:, j)
417  end do
418 
419  else if( flag == updatelag ) then
420 
421  call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv)
422  b(1:6, 1:nn*ndof) = 0.0d0
423  do j = 1, nn
424  b4 = ( bbar2(j, 1)-gderiv(j, 1) )/3.0d0
425  b6 = ( bbar2(j, 2)-gderiv(j, 2) )/3.0d0
426  b8 = ( bbar2(j, 3)-gderiv(j, 3) )/3.0d0
427  b(1, 3*j-2) = gderiv(j, 1)+b4
428  b(1, 3*j-1) = b6
429  b(1, 3*j ) = b8
430 
431  b(2, 3*j-2) = b4
432  b(2, 3*j-1) = gderiv(j, 2)+b6
433  b(2, 3*j ) = b8
434 
435  b(3, 3*j-2) = b4
436  b(3, 3*j-1) = b6
437  b(3, 3*j ) = gderiv(j, 3)+b8
438 
439  b(4, 3*j-2) = gderiv(j, 2)
440  b(4, 3*j-1) = gderiv(j, 1)
441  b(5, 3*j-1) = gderiv(j, 3)
442  b(5, 3*j ) = gderiv(j, 2)
443  b(6, 3*j-2) = gderiv(j, 3)
444  b(6, 3*j ) = gderiv(j, 1)
445  end do
446  end if
447 
448  !! calculate the Internal Force
449  wg = getweight(etype, lx)*det
450  qf(1:nn*ndof) = qf(1:nn*ndof) &
451  +matmul( gausses(lx)%stress(1:6), b(1:6, 1:nn*ndof) )*wg
452 
453  end do
454 
455  end subroutine update_c3d8bbar
456 
458  !----------------------------------------------------------------------*
459  subroutine tload_c3d8bbar &
460  (etype, nn, xx, yy, zz, tt, t0, &
461  gausses, vect, cdsys_id, coords)
462  !----------------------------------------------------------------------*
463 
464  use m_fstr
465  use mmechgauss
466  use m_matmatrix
467  use m_utilities
468 
469  !---------------------------------------------------------------------
470 
471  integer(kind=kint), parameter :: ndof = 3
472  integer(kind=kint), intent(in) :: etype, nn
473  type(tgaussstatus), intent(in) :: gausses(:)
474  real(kind=kreal), intent(in) :: xx(nn), yy(nn), zz(nn)
475  real(kind=kreal), intent(in) :: tt(nn), t0(nn)
476  real(kind=kreal), intent(out) :: vect(nn*ndof)
477  integer(kind=kint), intent(in) :: cdsys_ID
478  real(kind=kreal), intent(inout) :: coords(3, 3)
479 
480  !---------------------------------------------------------------------
481 
482  real(kind=kreal) :: alp, alp0, d(6, 6), b(6, ndof*nn)
483  real(kind=kreal) :: b4, b6, b8, det, ecoord(3, nn)
484  integer(kind=kint) :: j, LX, serr
485  real(kind=kreal) :: estrain(6), sgm(6), h(nn)
486  real(kind=kreal) :: naturalcoord(3), gderiv(nn, 3)
487  real(kind=kreal) :: wg, outa(1), ina(1), bbar(nn, 3), alpo(3), alpo0(3), coordsys(3, 3)
488  real(kind=kreal) :: tempc, temp0, tempc0, temp00, thermal_eps, tm(6,6)
489  logical :: ierr, matlaniso
490 
491  !---------------------------------------------------------------------
492 
493  matlaniso = .false.
494 
495  if( cdsys_id > 0 ) then ! cannot define aniso exapansion when no local coord defined
496  ina = tt(1)
497  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
498  if( .not. ierr ) matlaniso = .true.
499  end if
500 
501  vect(:) = 0.0d0
502 
503  ecoord(1, :) = xx(:)
504  ecoord(2, :) = yy(:)
505  ecoord(3, :) = zz(:)
506 
507  naturalcoord = 0.0d0
508  call getglobalderiv(etype, nn, naturalcoord, ecoord, det, bbar)
509  call getshapefunc( etype, naturalcoord, h(1:nn) )
510  tempc0 = dot_product( h(1:nn), tt(1:nn) )
511  temp00 = dot_product( h(1:nn), t0(1:nn) )
512 
513  ! LOOP FOR INTEGRATION POINTS
514  do lx = 1, numofquadpoints(etype)
515 
516  call getquadpoint( etype, lx, naturalcoord(:) )
517  call getshapefunc( etype, naturalcoord, h(1:nn) )
518  call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
519 
520  if( matlaniso ) then
521  call set_localcoordsys(coords, g_localcoordsys(cdsys_id), coordsys, serr)
522  if( serr == -1 ) stop "Fail to setup local coordinate"
523  if( serr == -2 ) then
524  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
525  end if
526  end if
527 
528  ! WEIGHT VALUE AT GAUSSIAN POINT
529  wg = getweight(etype, lx)*det
530  b(1:6, 1:nn*ndof) = 0.0d0
531  do j = 1, nn
532  b4 = ( bbar(j, 1)-gderiv(j, 1) )/3.0d0
533  b6 = ( bbar(j, 2)-gderiv(j, 2) )/3.0d0
534  b8 = ( bbar(j, 3)-gderiv(j, 3) )/3.0d0
535  b(1, 3*j-2) = gderiv(j, 1)+b4
536  b(1, 3*j-1) = b6
537  b(1, 3*j ) = b8
538 
539  b(2, 3*j-2) = b4
540  b(2, 3*j-1) = gderiv(j, 2)+b6
541  b(2, 3*j ) = b8
542 
543  b(3, 3*j-2) = b4
544  b(3, 3*j-1) = b6
545  b(3, 3*j ) = gderiv(j, 3)+b8
546 
547  b(4, 3*j-2) = gderiv(j, 2)
548  b(4, 3*j-1) = gderiv(j, 1)
549  b(5, 3*j-1) = gderiv(j, 3)
550  b(5, 3*j ) = gderiv(j, 2)
551  b(6, 3*j-2) = gderiv(j, 3)
552  b(6, 3*j ) = gderiv(j, 1)
553  end do
554 
555  tempc = dot_product( h(1:nn),tt(1:nn) )
556  temp0 = dot_product( h(1:nn),t0(1:nn) )
557 
558  call matlmatrix( gausses(lx), d3, d, 1.d0, 0.0d0, coordsys, tempc )
559 
560  ina(1) = tempc
561  if( matlaniso ) then
562  call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo(:), ierr, ina )
563  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
564  else
565  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
566  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
567  alp = outa(1)
568  end if
569  ina(1) = temp0
570  if( matlaniso ) then
571  call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo0(:), ierr, ina )
572  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
573  else
574  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
575  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
576  alp0 = outa(1)
577  end if
578 
579  !**
580  !** THERMAL strain
581  !**
582  if( matlaniso ) then
583  do j = 1, 3
584  estrain(j) = alpo(j)*(tempc0-ref_temp)-alpo0(j)*(temp00-ref_temp)
585  end do
586  estrain(4:6) = 0.0d0
587  call transformation(coordsys, tm)
588  estrain(:) = matmul( estrain(:), tm ) ! to global coord
589  estrain(4:6) = estrain(4:6)*2.0d0
590  else
591  thermal_eps = alp*(tempc0-ref_temp)-alp0*(temp00-ref_temp)
592  estrain(1:3) = thermal_eps
593  estrain(4:6) = 0.0d0
594  end if
595 
596  !**
597  !** SET SGM {s}=[D]{e}
598  !**
599  sgm(:) = matmul( d(:, :), estrain(:) )
600 
601  !**
602  !** CALCULATE LOAD {F}=[B]T{e}
603  !**
604  vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(:), b(:, :) )*wg
605 
606  end do
607 
608  end subroutine tload_c3d8bbar
609 
610 
611 end module m_static_lib_c3d8
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
Definition: element.f90:640
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:482
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:734
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:528
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:445
Definition: hecmw.f90:6
This modules defines common structures for fem analysis.
type(tlocalcoordsys), dimension(:), pointer, save g_localcoordsys
subroutine set_localcoordsys(coords, coordsys, outsys, ierr)
setup of coordinate system
This module provide functions for elastoplastic calculation.
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.f90:120
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
subroutine matlmatrix(gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp)
Calculate constituive matrix.
This module provide common functions of Solid elements.
subroutine geomat_c3(stress, mat)
subroutine cal_thermal_expansion_c3(tt0, ttc, material, coordsys, matlaniso, EPSTH)
subroutine update_stress3d(flag, gauss, rot, dstrain, F, coordsys, time, tincr, ttc, tt0, ttn)
This module contains several strategy to free locking problem in Eight-node hexagonal element.
subroutine tload_c3d8bbar(etype, nn, XX, YY, ZZ, TT, T0, gausses, VECT, cdsys_ID, coords)
This subroutien calculate thermal loading.
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.
subroutine stf_c3d8bbar(etype, nn, ecoord, gausses, stiff, cdsys_ID, coords, time, tincr, u, temperature)
This subroutine calculate stiff matrix using b-bar method.
This module provides aux functions.
Definition: utilities.f90:6
subroutine transformation(jacob, tm)
Definition: utilities.f90:339
This module provides functions for hyperelastic calculation.
Definition: Hyperelastic.f90:6
This module summarizes all infomation of material properties.
Definition: material.f90:6
integer(kind=kint), parameter totallag
Definition: material.f90:14
character(len=dict_key_length) mc_orthoexp
Definition: material.f90:123
integer(kind=kint), parameter infinitesimal
Definition: material.f90:13
integer(kind=kint), parameter updatelag
Definition: material.f90:15
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13