FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
static_LIB_Fbar.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  real(kind=kreal), private, parameter :: i33(3,3) = &
15  & reshape( (/1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 1.d0/), (/3,3/) )
16 
17 contains
18 
19 
21  !----------------------------------------------------------------------*
22  subroutine stf_c3d8fbar &
23  (etype, nn, ecoord, gausses, stiff, cdsys_id, coords, &
24  time, tincr, u, temperature)
25  !----------------------------------------------------------------------*
26 
27  use mmechgauss
28  use m_matmatrix
29  use m_common_struct
30  use m_static_lib_3d, only: geomat_c3
31  use m_utilities
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, p, q, lx, serr
55  real(kind=kreal) :: naturalcoord(3)
56  real(kind=kreal) :: gdispderiv(3, 3)
57  real(kind=kreal) :: b1(6, ndof*nn)
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) :: coordsys(3, 3)
61 
62  real(kind=kreal) :: elem0(3,nn), elem1(3, nn), gderiv1(nn,ndof), b2(6, ndof*nn), z1(3,2)
63  real(kind=kreal) :: v0, jacob, jacob_ave, gderiv1_ave(nn,ndof)
64  real(kind=kreal) :: gderiv2_ave(ndof*nn,ndof*nn)
65  real(kind=kreal) :: fbar(3,3), jratio(8), coeff, sff, dstrain(6), ddfs, fs(3,3), gfs(3,2)
66 
67  !---------------------------------------------------------------------
68 
69  stiff(:, :) = 0.0d0
70  ! we suppose the same material type in the element
71  flag = gausses(1)%pMaterial%nlgeom_flag
72  if( .not. present(u) ) flag = infinitesimal ! enforce to infinitesimal deformation analysis
73  elem(:, :) = ecoord(:, :)
74  elem0(:, :) = ecoord(:, :)
75  if( flag == updatelag ) elem(:, :) = ecoord(:, :)+u(:, :)
76  elem1(:, :) = ecoord(:, :)+u(:, :)
77 
78  !cal volumetric average of J=detF and dN/dx
79  v0 = 0.d0
80  jacob_ave = 0.d0
81  gderiv1_ave(1:nn,1:ndof) = 0.d0
82  gderiv2_ave(1:ndof*nn,1:ndof*nn) = 0.d0
83  do lx = 1, numofquadpoints(etype)
84  call getquadpoint( etype, lx, naturalcoord(:) )
85  call getglobalderiv( etype, nn, naturalcoord, elem0, det, gderiv)
86  wg = getweight(etype, lx)*det
87  if( flag == infinitesimal ) then
88  jacob = 1.d0
89  gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + jacob*wg*gderiv(1:nn, 1:ndof)
90  else
91  gdispderiv(1:3, 1:3) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
92  jacob = determinant33( i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof) )
93  jratio(lx) = dsign(dabs(jacob)**(-1.d0/3.d0),jacob)
94  call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv1)
95  gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + jacob*wg*gderiv1(1:nn, 1:ndof)
96  do p=1,nn
97  do q=1,nn
98  do i=1,ndof
99  do j=1,ndof
100  gderiv2_ave(3*p-3+i,3*q-3+j) = gderiv2_ave(3*p-3+i,3*q-3+j) &
101  & + jacob*wg*(gderiv1(p,i)*gderiv1(q,j)-gderiv1(q,i)*gderiv1(p,j))
102  end do
103  end do
104  end do
105  end do
106  endif
107  v0 = v0 + wg
108  jacob_ave = jacob_ave + jacob*wg
109  enddo
110  if( dabs(v0) > 1.d-12 ) then
111  if( dabs(jacob_ave) < 1.d-12 ) stop 'Error in Update_C3D8Fbar: too small average jacob'
112  jacob_ave = jacob_ave/v0
113  jratio(1:8) = (dsign(dabs(jacob_ave)**(1.d0/3.d0),jacob_ave))*jratio(1:8) !Jratio(1:8) = (jacob_ave**(1.d0/3.d0))*Jratio(1:8)
114  gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
115  gderiv2_ave(1:ndof*nn,1:ndof*nn) = gderiv2_ave(1:ndof*nn,1:ndof*nn)/(v0*jacob_ave)
116  else
117  stop 'Error in Update_C3D8Fbar: too small element volume'
118  end if
119 
120  do lx = 1, numofquadpoints(etype)
121 
122  call getquadpoint( etype, lx, naturalcoord(:) )
123  call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
124 
125  if( cdsys_id > 0 ) then
126  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
127  if( serr == -1 ) stop "Fail to setup local coordinate"
128  if( serr == -2 ) then
129  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
130  end if
131  end if
132 
133  if( present(temperature) ) then
134  call getshapefunc( etype, naturalcoord, spfunc )
135  temp = dot_product( temperature, spfunc )
136  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
137  else
138  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys )
139  end if
140 
141  if( flag == updatelag ) then
142  call geomat_c3( gausses(lx)%stress, mat )
143  d(:, :) = d(:, :)-mat
144  endif
145 
146  wg = getweight(etype, lx)*det
147  b(1:6, 1:nn*ndof) = 0.0d0
148  do j = 1, nn
149  b(1, 3*j-2) = gderiv(j, 1)
150  b(2, 3*j-1) = gderiv(j, 2)
151  b(3, 3*j ) = gderiv(j, 3)
152  b(4, 3*j-2) = gderiv(j, 2)
153  b(4, 3*j-1) = gderiv(j, 1)
154  b(5, 3*j-1) = gderiv(j, 3)
155  b(5, 3*j ) = gderiv(j, 2)
156  b(6, 3*j-2) = gderiv(j, 3)
157  b(6, 3*j ) = gderiv(j, 1)
158  end do
159 
160  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
161  if( flag == infinitesimal ) then
162  b2(1:6, 1:nn*ndof) = 0.0d0
163  do j = 1, nn
164  z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
165  b2(1,3*j-2:3*j) = z1(1:3,1)
166  b2(2,3*j-2:3*j) = z1(1:3,1)
167  b2(3,3*j-2:3*j) = z1(1:3,1)
168  end do
169 
170  ! BL = Jratio*(BL0 + BL1)+BL2
171  do j = 1, nn*ndof
172  b(1:3, j) = b(1:3,j)+b2(1:3,j)
173  end do
174 
175  else if( flag == totallag ) then
176  gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
177  fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
178  call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv1)
179 
180  ! ---dudx(i,j) ==> gdispderiv(i,j)
181  b1(1:6, 1:nn*ndof) = 0.0d0
182  do j = 1, nn
183  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
184  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
185  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
186  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
187  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
188  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
189  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
190  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
191  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
192  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
193  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
194  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
195  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
196  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
197  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
198  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
199  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
200  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
201  end do
202 
203  dstrain(1) = 0.5d0*(dot_product(fbar(1:3,1),fbar(1:3,1))-1.d0)
204  dstrain(2) = 0.5d0*(dot_product(fbar(1:3,2),fbar(1:3,2))-1.d0)
205  dstrain(3) = 0.5d0*(dot_product(fbar(1:3,3),fbar(1:3,3))-1.d0)
206  dstrain(4) = dot_product(fbar(1:3,1),fbar(1:3,2))
207  dstrain(5) = dot_product(fbar(1:3,2),fbar(1:3,3))
208  dstrain(6) = dot_product(fbar(1:3,3),fbar(1:3,1))
209 
210  b2(1:6, 1:nn*ndof) = 0.0d0
211  do j = 1, nn
212  z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
213  b2(1,3*j-2:3*j) = (2.d0*dstrain(1)+1.d0)*z1(1:3,1)
214  b2(2,3*j-2:3*j) = (2.d0*dstrain(2)+1.d0)*z1(1:3,1)
215  b2(3,3*j-2:3*j) = (2.d0*dstrain(3)+1.d0)*z1(1:3,1)
216  b2(4,3*j-2:3*j) = 2.d0*dstrain(4)*z1(1:3,1)
217  b2(5,3*j-2:3*j) = 2.d0*dstrain(5)*z1(1:3,1)
218  b2(6,3*j-2:3*j) = 2.d0*dstrain(6)*z1(1:3,1)
219  end do
220 
221  ! BL = Jratio*(BL0 + BL1)+BL2
222  do j = 1, nn*ndof
223  b(:, j) = jratio(lx)*jratio(lx)*(b(:,j)+b1(:,j))+b2(:,j)
224  end do
225 
226  else if( flag == updatelag ) then
227  wg = (jratio(lx)**3.d0)*getweight(etype, lx)*det
228 
229  b2(1:3, 1:nn*ndof) = 0.0d0
230  do j = 1, nn
231  z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
232  b2(1, 3*j-2:3*j) = z1(1:3,1)
233  b2(2, 3*j-2:3*j) = z1(1:3,1)
234  b2(3, 3*j-2:3*j) = z1(1:3,1)
235  end do
236 
237  do j = 1, nn*ndof
238  b(1:3, j) = b(1:3,j)+b2(1:3,j)
239  end do
240 
241  end if
242 
243  db(1:6, 1:nn*ndof) = matmul( d, b(1:6, 1:nn*ndof) )
244  forall( i=1:nn*ndof, j=1:nn*ndof )
245  stiff(i, j) = stiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
246  end forall
247 
248  ! calculate the initial stress matrix(1): dFbar*dFbar*Stress
249  if( flag == totallag .or. flag == updatelag ) then
250  stress(1:6) = gausses(lx)%stress
251  if( flag == totallag ) then
252  coeff = jratio(lx)
253  sff = dot_product(stress(1:6),dstrain(1:6))
254  else if( flag == updatelag ) then
255  coeff = 1.d0
256  sff = stress(1)+stress(2)+stress(3)
257  gderiv1 = gderiv
258  fbar(1:3,1:3) = i33(1:3,1:3)
259  end if
260 
261  bn(1:9, 1:nn*ndof) = 0.0d0
262  do j = 1, nn
263  bn(1, 3*j-2) = coeff*gderiv(j, 1)
264  bn(2, 3*j-1) = coeff*gderiv(j, 1)
265  bn(3, 3*j ) = coeff*gderiv(j, 1)
266  bn(4, 3*j-2) = coeff*gderiv(j, 2)
267  bn(5, 3*j-1) = coeff*gderiv(j, 2)
268  bn(6, 3*j ) = coeff*gderiv(j, 2)
269  bn(7, 3*j-2) = coeff*gderiv(j, 3)
270  bn(8, 3*j-1) = coeff*gderiv(j, 3)
271  bn(9, 3*j ) = coeff*gderiv(j, 3)
272  z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
273  bn(1, 3*j-2:3*j) = bn(1, 3*j-2:3*j) + fbar(1,1)*z1(1:3,1)
274  bn(2, 3*j-2:3*j) = bn(2, 3*j-2:3*j) + fbar(2,1)*z1(1:3,1)
275  bn(3, 3*j-2:3*j) = bn(3, 3*j-2:3*j) + fbar(3,1)*z1(1:3,1)
276  bn(4, 3*j-2:3*j) = bn(4, 3*j-2:3*j) + fbar(1,2)*z1(1:3,1)
277  bn(5, 3*j-2:3*j) = bn(5, 3*j-2:3*j) + fbar(2,2)*z1(1:3,1)
278  bn(6, 3*j-2:3*j) = bn(6, 3*j-2:3*j) + fbar(3,2)*z1(1:3,1)
279  bn(7, 3*j-2:3*j) = bn(7, 3*j-2:3*j) + fbar(1,3)*z1(1:3,1)
280  bn(8, 3*j-2:3*j) = bn(8, 3*j-2:3*j) + fbar(2,3)*z1(1:3,1)
281  bn(9, 3*j-2:3*j) = bn(9, 3*j-2:3*j) + fbar(3,3)*z1(1:3,1)
282  end do
283  smat(:, :) = 0.0d0
284  do j= 1, 3
285  smat(j , j ) = stress(1)
286  smat(j , j+3) = stress(4)
287  smat(j , j+6) = stress(6)
288  smat(j+3, j ) = stress(4)
289  smat(j+3, j+3) = stress(2)
290  smat(j+3, j+6) = stress(5)
291  smat(j+6, j ) = stress(6)
292  smat(j+6, j+3) = stress(5)
293  smat(j+6, j+6) = stress(3)
294  end do
295  sbn(1:9, 1:nn*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:nn*ndof) )
296  forall( i=1:nn*ndof, j=1:nn*ndof )
297  stiff(i, j) = stiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
298  end forall
299 
300  ! calculate the initial stress matrix(2): d(dFbar)*Stress
301  fs(1,1) = fbar(1,1)*stress(1)+fbar(1,2)*stress(4)+fbar(1,3)*stress(6)
302  fs(1,2) = fbar(1,1)*stress(4)+fbar(1,2)*stress(2)+fbar(1,3)*stress(5)
303  fs(1,3) = fbar(1,1)*stress(6)+fbar(1,2)*stress(5)+fbar(1,3)*stress(3)
304  fs(2,1) = fbar(2,1)*stress(1)+fbar(2,2)*stress(4)+fbar(2,3)*stress(6)
305  fs(2,2) = fbar(2,1)*stress(4)+fbar(2,2)*stress(2)+fbar(2,3)*stress(5)
306  fs(2,3) = fbar(2,1)*stress(6)+fbar(2,2)*stress(5)+fbar(2,3)*stress(3)
307  fs(3,1) = fbar(3,1)*stress(1)+fbar(3,2)*stress(4)+fbar(3,3)*stress(6)
308  fs(3,2) = fbar(3,1)*stress(4)+fbar(3,2)*stress(2)+fbar(3,3)*stress(5)
309  fs(3,3) = fbar(3,1)*stress(6)+fbar(3,2)*stress(5)+fbar(3,3)*stress(3)
310  do i=1,nn
311  z1(1:3,1) = (gderiv1_ave(i,1:3)-gderiv1(i,1:3))/3.d0
312  gfs(1:3,1) = coeff*matmul(fs,gderiv(i,1:3))
313  do j=1,nn
314  z1(1:3,2) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
315  gfs(1:3,2) = coeff*matmul(fs,gderiv(j,1:3))
316  do p=1,ndof
317  do q=1,ndof
318  ddfs = z1(p,1)*z1(q,2)
319  ddfs = ddfs + (gderiv2_ave(3*i-3+p,3*j-3+q)-gderiv1_ave(i,p)*gderiv1_ave(j,q))/3.d0
320  ddfs = ddfs + gderiv1(i,q)*gderiv1(j,p)/3.d0
321  ddfs = sff*ddfs + z1(p,1)*gfs(q,2)+z1(q,2)*gfs(p,1)
322  stiff(3*i-3+p, 3*j-3+q) = stiff(3*i-3+p, 3*j-3+q) + ddfs*wg
323  end do
324  end do
325  end do
326  end do
327  end if
328 
329 
330  end do ! gauss roop
331 
332  end subroutine stf_c3d8fbar
333 
334 
336  !----------------------------------------------------------------------*
337  subroutine update_c3d8fbar &
338  (etype, nn, ecoord, u, du, cdsys_id, coords, &
339  qf ,gausses, iter, time, tincr, tt,t0, tn )
340  !----------------------------------------------------------------------*
341 
342  use m_fstr
343  use mmaterial
344  use mmechgauss
345  use m_matmatrix
346  use m_elastoplastic
347  use mhyperelastic
348  use m_utilities
349  use m_static_lib_3d
350 
351  !---------------------------------------------------------------------
352 
353  integer(kind=kint), intent(in) :: etype
354  integer(kind=kint), intent(in) :: nn
355  real(kind=kreal), intent(in) :: ecoord(3, nn)
356  real(kind=kreal), intent(in) :: u(3, nn)
357  real(kind=kreal), intent(in) :: du(3, nn)
358  integer(kind=kint), intent(in) :: cdsys_id
359  real(kind=kreal), intent(inout) :: coords(3, 3)
360  real(kind=kreal), intent(out) :: qf(nn*3)
361  type(tgaussstatus), intent(inout) :: gausses(:)
362  integer, intent(in) :: iter
363  real(kind=kreal), intent(in) :: time
364  real(kind=kreal), intent(in) :: tincr
365  real(kind=kreal), intent(in), optional :: tt(nn)
366  real(kind=kreal), intent(in), optional :: t0(nn)
367  real(kind=kreal), intent(in), optional :: tn(nn)
368 
369  !---------------------------------------------------------------------
370 
371  integer(kind=kint) :: flag
372  integer(kind=kint), parameter :: ndof = 3
373  real(kind=kreal) :: b(6, ndof*nn), b1(6, ndof*nn)
374  real(kind=kreal) :: gderiv(nn, 3), gdispderiv(3, 3), det, wg
375  integer(kind=kint) :: j, lx, serr
376  real(kind=kreal) :: naturalcoord(3), rot(3, 3), spfunc(nn)
377  real(kind=kreal) :: totaldisp(3, nn), elem(3, nn), elem1(3, nn), coordsys(3, 3)
378  real(kind=kreal) :: dstrain(6)
379  real(kind=kreal) :: dvol
380  real(kind=kreal) :: ttc, tt0, ttn, alpo(3), ina(1), epsth(6)
381  logical :: ierr, matlaniso
382 
383  real(kind=kreal) :: elem0(3,nn), gderiv1(nn,ndof), b2(6, ndof*nn), z1(3)
384  real(kind=kreal) :: v0, jacob, jacob_ave, gderiv1_ave(nn,ndof)
385  real(kind=kreal) :: fbar(3,3), jratio(8)
386  real(kind=kreal) :: jacob_ave05, gderiv05_ave(nn,ndof)
387 
388  !---------------------------------------------------------------------
389 
390  qf(:) = 0.0d0
391  ! we suppose the same material type in the element
392  flag = gausses(1)%pMaterial%nlgeom_flag
393  elem0(1:ndof,1:nn) = ecoord(1:ndof,1:nn)
394  totaldisp(:, :) = u(:, :)+du(:, :)
395  if( flag == infinitesimal ) then
396  elem(:, :) = ecoord(:, :)
397  elem1(:, :) = ecoord(:, :)
398  else if( flag == totallag ) then
399  elem(:, :) = ecoord(:, :)
400  elem1(:, :) = ( du(:, :)+u(:, :) )+ecoord(:, :)
401  else if( flag == updatelag ) then
402  elem(:, :) = ( 0.5d0*du(:, :)+u(:, :) )+ecoord(:, :)
403  elem1(:, :) = ( du(:, :)+u(:, :) )+ecoord(:, :)
404  totaldisp(:, :) = du(:, :)
405  end if
406 
407  matlaniso = .false.
408  if( cdsys_id > 0 .AND. present(tt) ) then
409  ina = tt(1)
410  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
411  if( .not. ierr ) matlaniso = .true.
412  end if
413 
414  !cal volumetric average of J=detF and dN/dx
415  v0 = 0.d0
416  jacob_ave = 0.d0
417  gderiv1_ave(1:nn,1:ndof) = 0.d0
418  if( flag == updatelag ) then
419  jacob_ave05 = 0.d0
420  gderiv05_ave(1:nn,1:ndof) = 0.d0
421  endif
422  do lx = 1, numofquadpoints(etype)
423  call getquadpoint( etype, lx, naturalcoord(:) )
424  call getglobalderiv( etype, nn, naturalcoord, elem0, det, gderiv)
425  wg = getweight(etype, lx)*det
426  if( flag == infinitesimal ) then
427  jacob = 1.d0
428  gderiv1(1:nn, 1:ndof) = gderiv(1:nn, 1:ndof)
429  else
430  gdispderiv(1:3, 1:3) = matmul( du(1:ndof, 1:nn)+u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
431  jacob = determinant33( i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof) )
432  jratio(lx) = dsign(dabs(jacob)**(-1.d0/3.d0),jacob) ! Jratio(LX) = jacob**(-1.d0/3.d0)
433 
434  call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv1)
435  endif
436  v0 = v0 + wg
437  jacob_ave = jacob_ave + jacob*wg
438  gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + jacob*wg*gderiv1(1:nn, 1:ndof)
439  if( flag == updatelag ) then
440  call getglobalderiv( etype, nn, naturalcoord, elem, det, gderiv1)
441  wg = getweight(etype, lx)*det
442  jacob_ave05 = jacob_ave05 + wg
443  gderiv05_ave(1:nn,1:ndof) = gderiv05_ave(1:nn,1:ndof) + wg*gderiv1(1:nn, 1:ndof)
444  endif
445  enddo
446  if( dabs(v0) > 1.d-12 ) then
447  if( dabs(jacob_ave) < 1.d-12 ) stop 'Error in Update_C3D8Fbar: too small average jacob'
448  jacob_ave = jacob_ave/v0
449  jratio(1:8) = (dsign(dabs(jacob_ave)**(1.d0/3.d0),jacob_ave))*jratio(1:8) !Jratio(1:8) = (jacob_ave**(1.d0/3.d0))*Jratio(1:8)
450  gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
451  if( flag == updatelag ) gderiv05_ave(1:nn,1:ndof) = gderiv05_ave(1:nn,1:ndof)/jacob_ave05
452  else
453  stop 'Error in Update_C3D8Fbar: too small element volume'
454  end if
455 
456  do lx = 1, numofquadpoints(etype)
457 
458  call getquadpoint( etype, lx, naturalcoord(:) )
459  call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
460 
461  if( cdsys_id > 0 ) then
462  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
463  if( serr == -1 ) stop "Fail to setup local coordinate"
464  if( serr == -2 ) then
465  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
466  end if
467  end if
468 
469  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
470 
471  ! ========================================================
472  ! UPDATE STRAIN and STRESS
473  ! ========================================================
474 
475  ! Thermal Strain
476  epsth = 0.0d0
477  if( present(tt) .AND. present(t0) ) then
478  call getshapefunc(etype, naturalcoord, spfunc)
479  ttc = dot_product(tt, spfunc)
480  tt0 = dot_product(t0, spfunc)
481  ttn = dot_product(tn, spfunc)
482  call cal_thermal_expansion_c3( tt0, ttc, gausses(lx)%pMaterial, coordsys, matlaniso, epsth )
483  end if
484 
485  ! Update strain
486  if( flag == infinitesimal ) then
487  dvol = dot_product( totaldisp(1,1:nn), gderiv1_ave(1:nn,1) ) !du1/dx1
488  dvol = dvol + dot_product( totaldisp(2,1:nn), gderiv1_ave(1:nn,2) ) !du2/dx2
489  dvol = dvol + dot_product( totaldisp(3,1:nn), gderiv1_ave(1:nn,3) ) !du3/dx3
490  dvol = (dvol-(gdispderiv(1,1)+gdispderiv(2,2)+gdispderiv(3,3)))/3.d0
491  dstrain(1) = gdispderiv(1, 1) + dvol
492  dstrain(2) = gdispderiv(2, 2) + dvol
493  dstrain(3) = gdispderiv(3, 3) + dvol
494  dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
495  dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
496  dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
497  dstrain(:) = dstrain(:)-epsth(:)
498 
499  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
500 
501  else if( flag == totallag ) then
502  fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
503 
504  ! Green-Lagrange strain
505  dstrain(1) = 0.5d0*(dot_product(fbar(1:3,1),fbar(1:3,1))-1.d0)
506  dstrain(2) = 0.5d0*(dot_product(fbar(1:3,2),fbar(1:3,2))-1.d0)
507  dstrain(3) = 0.5d0*(dot_product(fbar(1:3,3),fbar(1:3,3))-1.d0)
508  dstrain(4) = dot_product(fbar(1:3,1),fbar(1:3,2))
509  dstrain(5) = dot_product(fbar(1:3,2),fbar(1:3,3))
510  dstrain(6) = dot_product(fbar(1:3,3),fbar(1:3,1))
511  dstrain(:) = dstrain(:)-epsth(:)
512 
513  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
514 
515  else if( flag == updatelag ) then
516  dvol = dot_product( totaldisp(1,1:nn), gderiv05_ave(1:nn,1) ) !du1/dx1
517  dvol = dvol + dot_product( totaldisp(2,1:nn), gderiv05_ave(1:nn,2) ) !du2/dx2
518  dvol = dvol + dot_product( totaldisp(3,1:nn), gderiv05_ave(1:nn,3) ) !du3/dx3
519  dvol = (dvol-(gdispderiv(1,1)+gdispderiv(2,2)+gdispderiv(3,3)))/3.d0
520  dstrain(1) = gdispderiv(1, 1) + dvol
521  dstrain(2) = gdispderiv(2, 2) + dvol
522  dstrain(3) = gdispderiv(3, 3) + dvol
523  dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
524  dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
525  dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
526  dstrain(:) = dstrain(:)-epsth(:)
527 
528  rot = 0.0d0
529  rot(1, 2)= 0.5d0*( gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
530  rot(2, 3)= 0.5d0*( gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
531  rot(1, 3)= 0.5d0*( gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
532 
533  gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
534 
535  call getglobalderiv( etype, nn, naturalcoord, elem0, det, gderiv1)
536  gdispderiv(1:ndof, 1:ndof) = matmul( du(1:ndof, 1:nn)+u(1:ndof, 1:nn), gderiv1(1:nn, 1:ndof) )
537  fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
538 
539  end if
540 
541  ! Update stress
542  if( present(tt) .AND. present(t0) ) then
543  call update_stress3d( flag, gausses(lx), rot, dstrain, fbar, coordsys, time, tincr, ttc, tt0, ttn )
544  else
545  call update_stress3d( flag, gausses(lx), rot, dstrain, fbar, coordsys, time, tincr )
546  end if
547 
548  ! ========================================================
549  ! calculate the internal force ( equivalent nodal force )
550  ! ========================================================
551  ! Small strain
552  b(1:6, 1:nn*ndof) = 0.0d0
553  do j = 1, nn
554  b(1,3*j-2) = gderiv(j, 1)
555  b(2,3*j-1) = gderiv(j, 2)
556  b(3,3*j ) = gderiv(j, 3)
557  b(4,3*j-2) = gderiv(j, 2)
558  b(4,3*j-1) = gderiv(j, 1)
559  b(5,3*j-1) = gderiv(j, 3)
560  b(5,3*j ) = gderiv(j, 2)
561  b(6,3*j-2) = gderiv(j, 3)
562  b(6,3*j ) = gderiv(j, 1)
563  end do
564 
565  wg=getweight( etype, lx )*det
566  if( flag == infinitesimal ) then
567  gderiv1(1:nn, 1:ndof) = gderiv(1:nn, 1:ndof)
568 
569  b2(1:6, 1:nn*ndof) = 0.0d0
570  do j = 1, nn
571  z1(1:3) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
572  b2(1,3*j-2:3*j) = z1(1:3)
573  b2(2,3*j-2:3*j) = z1(1:3)
574  b2(3,3*j-2:3*j) = z1(1:3)
575  end do
576 
577  ! BL = BL0 + BL2
578  do j = 1, nn*ndof
579  b(:, j) = b(:,j)+b2(:,j)
580  end do
581 
582  else if( flag == totallag ) then
583 
584  b1(1:6, 1:nn*ndof) = 0.0d0
585  do j = 1, nn
586  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
587  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
588  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
589  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
590  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
591  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
592  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
593  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
594  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
595  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
596  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
597  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
598  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
599  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
600  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
601  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
602  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
603  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
604  end do
605 
606  call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv1)
607 
608  b2(1:6, 1:nn*ndof) = 0.0d0
609  do j = 1, nn
610  z1(1:3) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
611  b2(1,3*j-2:3*j) = (2.d0*dstrain(1)+1.d0)*z1(1:3)
612  b2(2,3*j-2:3*j) = (2.d0*dstrain(2)+1.d0)*z1(1:3)
613  b2(3,3*j-2:3*j) = (2.d0*dstrain(3)+1.d0)*z1(1:3)
614  b2(4,3*j-2:3*j) = 2.d0*dstrain(4)*z1(1:3)
615  b2(5,3*j-2:3*j) = 2.d0*dstrain(5)*z1(1:3)
616  b2(6,3*j-2:3*j) = 2.d0*dstrain(6)*z1(1:3)
617  end do
618 
619  ! BL = R^2*(BL0 + BL1)+BL2
620  do j = 1, nn*ndof
621  b(:, j) = jratio(lx)*jratio(lx)*(b(:,j)+b1(:,j))+b2(:,j)
622  end do
623 
624  else if( flag == updatelag ) then
625 
626  call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv)
627  wg = (jratio(lx)**3.d0)*getweight(etype, lx)*det
628  b(1:6, 1:nn*ndof) = 0.0d0
629  do j = 1, nn
630  b(1, 3*j-2) = gderiv(j, 1)
631  b(2, 3*j-1) = gderiv(j, 2)
632  b(3, 3*j ) = gderiv(j, 3)
633  b(4, 3*j-2) = gderiv(j, 2)
634  b(4, 3*j-1) = gderiv(j, 1)
635  b(5, 3*j-1) = gderiv(j, 3)
636  b(5, 3*j ) = gderiv(j, 2)
637  b(6, 3*j-2) = gderiv(j, 3)
638  b(6, 3*j ) = gderiv(j, 1)
639  end do
640 
641  b2(1:3, 1:nn*ndof) = 0.0d0
642  do j = 1, nn
643  z1(1:3) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
644  b2(1, 3*j-2:3*j) = z1(1:3)
645  b2(2, 3*j-2:3*j) = z1(1:3)
646  b2(3, 3*j-2:3*j) = z1(1:3)
647  end do
648 
649  do j = 1, nn*ndof
650  b(1:3, j) = b(1:3,j)+b2(1:3,j)
651  end do
652 
653  end if
654 
655  qf(1:nn*ndof) &
656  = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
657 
658  end do
659 
660  end subroutine update_c3d8fbar
661 
663  !----------------------------------------------------------------------*
664  subroutine tload_c3d8fbar &
665  (etype, nn, xx, yy, zz, tt, t0, &
666  gausses, vect, cdsys_id, coords)
667  !----------------------------------------------------------------------*
668 
669  use m_fstr
670  use mmechgauss
671  use m_matmatrix
672  use m_utilities
673 
674  !---------------------------------------------------------------------
675 
676  integer(kind=kint), parameter :: ndof = 3
677  integer(kind=kint), intent(in) :: etype, nn
678  type(tgaussstatus), intent(in) :: gausses(:)
679  real(kind=kreal), intent(in) :: xx(nn), yy(nn), zz(nn)
680  real(kind=kreal), intent(in) :: tt(nn), t0(nn)
681  real(kind=kreal), intent(out) :: vect(nn*ndof)
682  integer(kind=kint), intent(in) :: cdsys_ID
683  real(kind=kreal), intent(inout) :: coords(3, 3)
684 
685  !---------------------------------------------------------------------
686 
687  real(kind=kreal) :: alp, alp0, d(6, 6), b(6, ndof*nn)
688  real(kind=kreal) :: z1(3), det, ecoord(3, nn)
689  integer(kind=kint) :: j, LX, serr
690  real(kind=kreal) :: estrain(6), sgm(6), h(nn)
691  real(kind=kreal) :: naturalcoord(3), gderiv(nn, 3)
692  real(kind=kreal) :: wg, outa(1), ina(1), gderiv1_ave(nn, 3), alpo(3), alpo0(3), coordsys(3, 3)
693  real(kind=kreal) :: tempc, temp0, v0, jacob_ave, thermal_eps, tm(6,6)
694  logical :: ierr, matlaniso
695 
696  !---------------------------------------------------------------------
697 
698  matlaniso = .false.
699 
700  if( cdsys_id > 0 ) then ! cannot define aniso exapansion when no local coord defined
701  ina = tt(1)
702  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
703  if( .not. ierr ) matlaniso = .true.
704  end if
705 
706  vect(:) = 0.0d0
707 
708  ecoord(1, :) = xx(:)
709  ecoord(2, :) = yy(:)
710  ecoord(3, :) = zz(:)
711 
712  !cal volumetric average of J=detF and dN/dx
713  v0 = 0.d0
714  jacob_ave = 0.d0
715  gderiv1_ave(1:nn,1:ndof) = 0.d0
716  do lx = 1, numofquadpoints(etype)
717  call getquadpoint( etype, lx, naturalcoord(:) )
718  call getglobalderiv( etype, nn, naturalcoord, ecoord, det, gderiv)
719  wg = getweight(etype, lx)*det
720  v0 = v0 + wg
721  jacob_ave = jacob_ave + wg
722  gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + wg*gderiv(1:nn, 1:ndof)
723  enddo
724  if( dabs(v0) > 1.d-12 ) then
725  if( dabs(jacob_ave) < 1.d-12 ) stop 'Error in TLOAD_C3D8Fbar: too small average jacob'
726  jacob_ave = jacob_ave/v0
727  gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
728  else
729  stop 'Error in TLOAD_C3D8Fbar: too small element volume'
730  end if
731 
732  ! LOOP FOR INTEGRATION POINTS
733  do lx = 1, numofquadpoints(etype)
734 
735  call getquadpoint( etype, lx, naturalcoord(:) )
736  call getshapefunc( etype, naturalcoord, h(1:nn) )
737  call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
738 
739  if( matlaniso ) then
740  call set_localcoordsys(coords, g_localcoordsys(cdsys_id), coordsys, serr)
741  if( serr == -1 ) stop "Fail to setup local coordinate"
742  if( serr == -2 ) then
743  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
744  end if
745  end if
746 
747  ! WEIGHT VALUE AT GAUSSIAN POINT
748  wg = getweight(etype, lx)*det
749  b(1:6, 1:nn*ndof) = 0.0d0
750  do j = 1, nn
751  b(1,3*j-2) = gderiv(j, 1)
752  b(2,3*j-1) = gderiv(j, 2)
753  b(3,3*j ) = gderiv(j, 3)
754  b(4,3*j-2) = gderiv(j, 2)
755  b(4,3*j-1) = gderiv(j, 1)
756  b(5,3*j-1) = gderiv(j, 3)
757  b(5,3*j ) = gderiv(j, 2)
758  b(6,3*j-2) = gderiv(j, 3)
759  b(6,3*j ) = gderiv(j, 1)
760  end do
761 
762  do j = 1, nn
763  z1(1:3) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
764  b(1,3*j-2:3*j) = b(1,3*j-2:3*j)+z1(1:3)
765  b(2,3*j-2:3*j) = b(2,3*j-2:3*j)+z1(1:3)
766  b(3,3*j-2:3*j) = b(3,3*j-2:3*j)+z1(1:3)
767  end do
768 
769  tempc = dot_product( h(1:nn),tt(1:nn) )
770  temp0 = dot_product( h(1:nn),t0(1:nn) )
771 
772  call matlmatrix( gausses(lx), d3, d, 1.d0, 0.0d0, coordsys, tempc )
773 
774  ina(1) = tempc
775  if( matlaniso ) then
776  call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo(:), ierr, ina )
777  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
778  else
779  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
780  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
781  alp = outa(1)
782  end if
783  ina(1) = temp0
784  if( matlaniso ) then
785  call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo0(:), ierr, ina )
786  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
787  else
788  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
789  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
790  alp0 = outa(1)
791  end if
792 
793  !**
794  !** THERMAL strain
795  !**
796  if( matlaniso ) then
797  do j = 1, 3
798  estrain(j) = alpo(j)*(tempc-ref_temp)-alpo0(j)*(temp0-ref_temp)
799  end do
800  estrain(4:6) = 0.0d0
801  call transformation(coordsys, tm)
802  estrain(:) = matmul( estrain(:), tm ) ! to global coord
803  estrain(4:6) = estrain(4:6)*2.0d0
804  else
805  thermal_eps = alp*(tempc-ref_temp)-alp0*(temp0-ref_temp)
806  estrain(1:3) = thermal_eps
807  estrain(4:6) = 0.0d0
808  end if
809 
810  !**
811  !** SET SGM {s}=[D]{e}
812  !**
813  sgm(:) = matmul( d(:, :), estrain(:) )
814 
815  !**
816  !** CALCULATE LOAD {F}=[B]T{e}
817  !**
818  vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(:), b(:, :) )*wg
819 
820  end do
821 
822  end subroutine tload_c3d8fbar
823 
824 
825 end module m_static_lib_fbar
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 update_c3d8fbar(etype, nn, ecoord, u, du, cdsys_ID, coords, qf, gausses, iter, time, tincr, TT, T0, TN)
Update Strain stress of this element.
subroutine stf_c3d8fbar(etype, nn, ecoord, gausses, stiff, cdsys_ID, coords, time, tincr, u, temperature)
This subroutine calculate stiff matrix using F-bar method.
subroutine tload_c3d8fbar(etype, nn, XX, YY, ZZ, TT, T0, gausses, VECT, cdsys_ID, coords)
This subroutien calculate thermal loading.
This module provides aux functions.
Definition: utilities.f90:6
real(kind=kreal) function determinant33(XJ)
Compute determinant for 3*3 matrix.
Definition: utilities.f90:233
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