FrontISTR  5.7.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 !-------------------------------------------------------------------------------
7 module m_static_lib_fbar
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) :: 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  call getshapefunc( etype, naturalcoord, spfunc )
134  temp = dot_product( temperature, spfunc )
135  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
136 
137  if( flag == updatelag ) then
138  call geomat_c3( gausses(lx)%stress, mat )
139  d(:, :) = d(:, :)-mat
140  endif
141 
142  wg = getweight(etype, lx)*det
143  b(1:6, 1:nn*ndof) = 0.0d0
144  do j = 1, nn
145  b(1, 3*j-2) = gderiv(j, 1)
146  b(2, 3*j-1) = gderiv(j, 2)
147  b(3, 3*j ) = gderiv(j, 3)
148  b(4, 3*j-2) = gderiv(j, 2)
149  b(4, 3*j-1) = gderiv(j, 1)
150  b(5, 3*j-1) = gderiv(j, 3)
151  b(5, 3*j ) = gderiv(j, 2)
152  b(6, 3*j-2) = gderiv(j, 3)
153  b(6, 3*j ) = gderiv(j, 1)
154  end do
155 
156  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
157  if( flag == infinitesimal ) then
158  b2(1:6, 1:nn*ndof) = 0.0d0
159  do j = 1, nn
160  z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
161  b2(1,3*j-2:3*j) = z1(1:3,1)
162  b2(2,3*j-2:3*j) = z1(1:3,1)
163  b2(3,3*j-2:3*j) = z1(1:3,1)
164  end do
165 
166  ! BL = Jratio*(BL0 + BL1)+BL2
167  do j = 1, nn*ndof
168  b(1:3, j) = b(1:3,j)+b2(1:3,j)
169  end do
170 
171  else if( flag == totallag ) then
172  gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
173  fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
174  call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv1)
175 
176  ! ---dudx(i,j) ==> gdispderiv(i,j)
177  b1(1:6, 1:nn*ndof) = 0.0d0
178  do j = 1, nn
179  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
180  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
181  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
182  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
183  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
184  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
185  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
186  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
187  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
188  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
189  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
190  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
191  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
192  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
193  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
194  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
195  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
196  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
197  end do
198 
199  dstrain(1) = 0.5d0*(dot_product(fbar(1:3,1),fbar(1:3,1))-1.d0)
200  dstrain(2) = 0.5d0*(dot_product(fbar(1:3,2),fbar(1:3,2))-1.d0)
201  dstrain(3) = 0.5d0*(dot_product(fbar(1:3,3),fbar(1:3,3))-1.d0)
202  dstrain(4) = dot_product(fbar(1:3,1),fbar(1:3,2))
203  dstrain(5) = dot_product(fbar(1:3,2),fbar(1:3,3))
204  dstrain(6) = dot_product(fbar(1:3,3),fbar(1:3,1))
205 
206  b2(1:6, 1:nn*ndof) = 0.0d0
207  do j = 1, nn
208  z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
209  b2(1,3*j-2:3*j) = (2.d0*dstrain(1)+1.d0)*z1(1:3,1)
210  b2(2,3*j-2:3*j) = (2.d0*dstrain(2)+1.d0)*z1(1:3,1)
211  b2(3,3*j-2:3*j) = (2.d0*dstrain(3)+1.d0)*z1(1:3,1)
212  b2(4,3*j-2:3*j) = 2.d0*dstrain(4)*z1(1:3,1)
213  b2(5,3*j-2:3*j) = 2.d0*dstrain(5)*z1(1:3,1)
214  b2(6,3*j-2:3*j) = 2.d0*dstrain(6)*z1(1:3,1)
215  end do
216 
217  ! BL = Jratio*(BL0 + BL1)+BL2
218  do j = 1, nn*ndof
219  b(:, j) = jratio(lx)*jratio(lx)*(b(:,j)+b1(:,j))+b2(:,j)
220  end do
221 
222  else if( flag == updatelag ) then
223  wg = (jratio(lx)**3.d0)*getweight(etype, lx)*det
224 
225  b2(1:3, 1:nn*ndof) = 0.0d0
226  do j = 1, nn
227  z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
228  b2(1, 3*j-2:3*j) = z1(1:3,1)
229  b2(2, 3*j-2:3*j) = z1(1:3,1)
230  b2(3, 3*j-2:3*j) = z1(1:3,1)
231  end do
232 
233  do j = 1, nn*ndof
234  b(1:3, j) = b(1:3,j)+b2(1:3,j)
235  end do
236 
237  end if
238 
239  db(1:6, 1:nn*ndof) = matmul( d, b(1:6, 1:nn*ndof) )
240  do j=1,nn*ndof
241  do i=1,nn*ndof
242  stiff(i, j) = stiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
243  end do
244  end do
245 
246  ! calculate the initial stress matrix(1): dFbar*dFbar*Stress
247  if( flag == totallag .or. flag == updatelag ) then
248  stress(1:6) = gausses(lx)%stress
249  if( flag == totallag ) then
250  coeff = jratio(lx)
251  sff = dot_product(stress(1:6),dstrain(1:6))
252  else if( flag == updatelag ) then
253  coeff = 1.d0
254  sff = stress(1)+stress(2)+stress(3)
255  gderiv1 = gderiv
256  fbar(1:3,1:3) = i33(1:3,1:3)
257  end if
258 
259  bn(1:9, 1:nn*ndof) = 0.0d0
260  do j = 1, nn
261  bn(1, 3*j-2) = coeff*gderiv(j, 1)
262  bn(2, 3*j-1) = coeff*gderiv(j, 1)
263  bn(3, 3*j ) = coeff*gderiv(j, 1)
264  bn(4, 3*j-2) = coeff*gderiv(j, 2)
265  bn(5, 3*j-1) = coeff*gderiv(j, 2)
266  bn(6, 3*j ) = coeff*gderiv(j, 2)
267  bn(7, 3*j-2) = coeff*gderiv(j, 3)
268  bn(8, 3*j-1) = coeff*gderiv(j, 3)
269  bn(9, 3*j ) = coeff*gderiv(j, 3)
270  z1(1:3,1) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
271  bn(1, 3*j-2:3*j) = bn(1, 3*j-2:3*j) + fbar(1,1)*z1(1:3,1)
272  bn(2, 3*j-2:3*j) = bn(2, 3*j-2:3*j) + fbar(2,1)*z1(1:3,1)
273  bn(3, 3*j-2:3*j) = bn(3, 3*j-2:3*j) + fbar(3,1)*z1(1:3,1)
274  bn(4, 3*j-2:3*j) = bn(4, 3*j-2:3*j) + fbar(1,2)*z1(1:3,1)
275  bn(5, 3*j-2:3*j) = bn(5, 3*j-2:3*j) + fbar(2,2)*z1(1:3,1)
276  bn(6, 3*j-2:3*j) = bn(6, 3*j-2:3*j) + fbar(3,2)*z1(1:3,1)
277  bn(7, 3*j-2:3*j) = bn(7, 3*j-2:3*j) + fbar(1,3)*z1(1:3,1)
278  bn(8, 3*j-2:3*j) = bn(8, 3*j-2:3*j) + fbar(2,3)*z1(1:3,1)
279  bn(9, 3*j-2:3*j) = bn(9, 3*j-2:3*j) + fbar(3,3)*z1(1:3,1)
280  end do
281  smat(:, :) = 0.0d0
282  do j= 1, 3
283  smat(j , j ) = stress(1)
284  smat(j , j+3) = stress(4)
285  smat(j , j+6) = stress(6)
286  smat(j+3, j ) = stress(4)
287  smat(j+3, j+3) = stress(2)
288  smat(j+3, j+6) = stress(5)
289  smat(j+6, j ) = stress(6)
290  smat(j+6, j+3) = stress(5)
291  smat(j+6, j+6) = stress(3)
292  end do
293  sbn(1:9, 1:nn*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:nn*ndof) )
294  do j=1,nn*ndof
295  do i=1,nn*ndof
296  stiff(i, j) = stiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
297  end do
298  end do
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) :: tt(nn)
366  real(kind=kreal), intent(in) :: t0(nn)
367  real(kind=kreal), intent(in) :: 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  ina = tt(1)
409  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
410  if( .not. ierr ) matlaniso = .true.
411 
412  !cal volumetric average of J=detF and dN/dx
413  v0 = 0.d0
414  jacob_ave = 0.d0
415  gderiv1_ave(1:nn,1:ndof) = 0.d0
416  if( flag == updatelag ) then
417  jacob_ave05 = 0.d0
418  gderiv05_ave(1:nn,1:ndof) = 0.d0
419  endif
420  do lx = 1, numofquadpoints(etype)
421  call getquadpoint( etype, lx, naturalcoord(:) )
422  call getglobalderiv( etype, nn, naturalcoord, elem0, det, gderiv)
423  wg = getweight(etype, lx)*det
424  if( flag == infinitesimal ) then
425  jacob = 1.d0
426  gderiv1(1:nn, 1:ndof) = gderiv(1:nn, 1:ndof)
427  else
428  gdispderiv(1:3, 1:3) = matmul( du(1:ndof, 1:nn)+u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
429  jacob = determinant33( i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof) )
430  jratio(lx) = dsign(dabs(jacob)**(-1.d0/3.d0),jacob) ! Jratio(LX) = jacob**(-1.d0/3.d0)
431 
432  call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv1)
433  endif
434  v0 = v0 + wg
435  jacob_ave = jacob_ave + jacob*wg
436  gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + jacob*wg*gderiv1(1:nn, 1:ndof)
437  if( flag == updatelag ) then
438  call getglobalderiv( etype, nn, naturalcoord, elem, det, gderiv1)
439  wg = getweight(etype, lx)*det
440  jacob_ave05 = jacob_ave05 + wg
441  gderiv05_ave(1:nn,1:ndof) = gderiv05_ave(1:nn,1:ndof) + wg*gderiv1(1:nn, 1:ndof)
442  endif
443  enddo
444  if( dabs(v0) > 1.d-12 ) then
445  if( dabs(jacob_ave) < 1.d-12 ) stop 'Error in Update_C3D8Fbar: too small average jacob'
446  jacob_ave = jacob_ave/v0
447  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)
448  gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
449  if( flag == updatelag ) gderiv05_ave(1:nn,1:ndof) = gderiv05_ave(1:nn,1:ndof)/jacob_ave05
450  else
451  stop 'Error in Update_C3D8Fbar: too small element volume'
452  end if
453 
454  do lx = 1, numofquadpoints(etype)
455 
456  call getquadpoint( etype, lx, naturalcoord(:) )
457  call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
458 
459  if( cdsys_id > 0 ) then
460  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
461  if( serr == -1 ) stop "Fail to setup local coordinate"
462  if( serr == -2 ) then
463  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
464  end if
465  end if
466 
467  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
468 
469  ! ========================================================
470  ! UPDATE STRAIN and STRESS
471  ! ========================================================
472 
473  ! Thermal Strain
474  call getshapefunc(etype, naturalcoord, spfunc)
475  ttc = dot_product(tt, spfunc)
476  tt0 = dot_product(t0, spfunc)
477  ttn = dot_product(tn, spfunc)
478  call cal_thermal_expansion_c3( tt0, ttc, gausses(lx)%pMaterial, coordsys, matlaniso, epsth )
479 
480  ! Update strain
481  if( flag == infinitesimal ) then
482  dvol = dot_product( totaldisp(1,1:nn), gderiv1_ave(1:nn,1) ) !du1/dx1
483  dvol = dvol + dot_product( totaldisp(2,1:nn), gderiv1_ave(1:nn,2) ) !du2/dx2
484  dvol = dvol + dot_product( totaldisp(3,1:nn), gderiv1_ave(1:nn,3) ) !du3/dx3
485  dvol = (dvol-(gdispderiv(1,1)+gdispderiv(2,2)+gdispderiv(3,3)))/3.d0
486  dstrain(1) = gdispderiv(1, 1) + dvol
487  dstrain(2) = gdispderiv(2, 2) + dvol
488  dstrain(3) = gdispderiv(3, 3) + dvol
489  dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
490  dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
491  dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
492  dstrain(:) = dstrain(:)-epsth(:)
493 
494  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
495 
496  else if( flag == totallag ) then
497  fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
498 
499  ! Green-Lagrange strain
500  dstrain(1) = 0.5d0*(dot_product(fbar(1:3,1),fbar(1:3,1))-1.d0)
501  dstrain(2) = 0.5d0*(dot_product(fbar(1:3,2),fbar(1:3,2))-1.d0)
502  dstrain(3) = 0.5d0*(dot_product(fbar(1:3,3),fbar(1:3,3))-1.d0)
503  dstrain(4) = dot_product(fbar(1:3,1),fbar(1:3,2))
504  dstrain(5) = dot_product(fbar(1:3,2),fbar(1:3,3))
505  dstrain(6) = dot_product(fbar(1:3,3),fbar(1:3,1))
506  dstrain(:) = dstrain(:)-epsth(:)
507 
508  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
509 
510  else if( flag == updatelag ) then
511  dvol = dot_product( totaldisp(1,1:nn), gderiv05_ave(1:nn,1) ) !du1/dx1
512  dvol = dvol + dot_product( totaldisp(2,1:nn), gderiv05_ave(1:nn,2) ) !du2/dx2
513  dvol = dvol + dot_product( totaldisp(3,1:nn), gderiv05_ave(1:nn,3) ) !du3/dx3
514  dvol = (dvol-(gdispderiv(1,1)+gdispderiv(2,2)+gdispderiv(3,3)))/3.d0
515  dstrain(1) = gdispderiv(1, 1) + dvol
516  dstrain(2) = gdispderiv(2, 2) + dvol
517  dstrain(3) = gdispderiv(3, 3) + dvol
518  dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
519  dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
520  dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
521  dstrain(:) = dstrain(:)-epsth(:)
522 
523  rot = 0.0d0
524  rot(1, 2)= 0.5d0*( gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
525  rot(2, 3)= 0.5d0*( gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
526  rot(1, 3)= 0.5d0*( gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
527 
528  gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
529 
530  call getglobalderiv( etype, nn, naturalcoord, elem0, det, gderiv1)
531  gdispderiv(1:ndof, 1:ndof) = matmul( du(1:ndof, 1:nn)+u(1:ndof, 1:nn), gderiv1(1:nn, 1:ndof) )
532  fbar(1:ndof, 1:ndof) = jratio(lx)*(i33(1:ndof,1:ndof) + gdispderiv(1:ndof, 1:ndof))
533 
534  end if
535 
536  ! Update stress
537  call update_stress3d( flag, gausses(lx), rot, dstrain, fbar, coordsys, time, tincr, ttc, tt0, ttn )
538 
539  ! ========================================================
540  ! calculate the internal force ( equivalent nodal force )
541  ! ========================================================
542  ! Small strain
543  b(1:6, 1:nn*ndof) = 0.0d0
544  do j = 1, nn
545  b(1,3*j-2) = gderiv(j, 1)
546  b(2,3*j-1) = gderiv(j, 2)
547  b(3,3*j ) = gderiv(j, 3)
548  b(4,3*j-2) = gderiv(j, 2)
549  b(4,3*j-1) = gderiv(j, 1)
550  b(5,3*j-1) = gderiv(j, 3)
551  b(5,3*j ) = gderiv(j, 2)
552  b(6,3*j-2) = gderiv(j, 3)
553  b(6,3*j ) = gderiv(j, 1)
554  end do
555 
556  wg=getweight( etype, lx )*det
557  if( flag == infinitesimal ) then
558  gderiv1(1:nn, 1:ndof) = gderiv(1:nn, 1:ndof)
559 
560  b2(1:6, 1:nn*ndof) = 0.0d0
561  do j = 1, nn
562  z1(1:3) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
563  b2(1,3*j-2:3*j) = z1(1:3)
564  b2(2,3*j-2:3*j) = z1(1:3)
565  b2(3,3*j-2:3*j) = z1(1:3)
566  end do
567 
568  ! BL = BL0 + BL2
569  do j = 1, nn*ndof
570  b(:, j) = b(:,j)+b2(:,j)
571  end do
572 
573  else if( flag == totallag ) then
574 
575  b1(1:6, 1:nn*ndof) = 0.0d0
576  do j = 1, nn
577  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
578  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
579  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
580  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
581  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
582  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
583  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
584  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
585  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
586  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
587  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
588  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
589  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
590  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
591  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
592  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
593  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
594  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
595  end do
596 
597  call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv1)
598 
599  b2(1:6, 1:nn*ndof) = 0.0d0
600  do j = 1, nn
601  z1(1:3) = (gderiv1_ave(j,1:3)-gderiv1(j,1:3))/3.d0
602  b2(1,3*j-2:3*j) = (2.d0*dstrain(1)+1.d0)*z1(1:3)
603  b2(2,3*j-2:3*j) = (2.d0*dstrain(2)+1.d0)*z1(1:3)
604  b2(3,3*j-2:3*j) = (2.d0*dstrain(3)+1.d0)*z1(1:3)
605  b2(4,3*j-2:3*j) = 2.d0*dstrain(4)*z1(1:3)
606  b2(5,3*j-2:3*j) = 2.d0*dstrain(5)*z1(1:3)
607  b2(6,3*j-2:3*j) = 2.d0*dstrain(6)*z1(1:3)
608  end do
609 
610  ! BL = R^2*(BL0 + BL1)+BL2
611  do j = 1, nn*ndof
612  b(:, j) = jratio(lx)*jratio(lx)*(b(:,j)+b1(:,j))+b2(:,j)
613  end do
614 
615  else if( flag == updatelag ) then
616 
617  call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv)
618  wg = (jratio(lx)**3.d0)*getweight(etype, lx)*det
619  b(1:6, 1:nn*ndof) = 0.0d0
620  do j = 1, nn
621  b(1, 3*j-2) = gderiv(j, 1)
622  b(2, 3*j-1) = gderiv(j, 2)
623  b(3, 3*j ) = gderiv(j, 3)
624  b(4, 3*j-2) = gderiv(j, 2)
625  b(4, 3*j-1) = gderiv(j, 1)
626  b(5, 3*j-1) = gderiv(j, 3)
627  b(5, 3*j ) = gderiv(j, 2)
628  b(6, 3*j-2) = gderiv(j, 3)
629  b(6, 3*j ) = gderiv(j, 1)
630  end do
631 
632  b2(1:3, 1:nn*ndof) = 0.0d0
633  do j = 1, nn
634  z1(1:3) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
635  b2(1, 3*j-2:3*j) = z1(1:3)
636  b2(2, 3*j-2:3*j) = z1(1:3)
637  b2(3, 3*j-2:3*j) = z1(1:3)
638  end do
639 
640  do j = 1, nn*ndof
641  b(1:3, j) = b(1:3,j)+b2(1:3,j)
642  end do
643 
644  end if
645 
646  qf(1:nn*ndof) &
647  = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
648 
649  ! integrate strain energy
650  gausses(lx)%strain_energy = gausses(lx)%strain_energy*wg
651  end do
652 
653  end subroutine update_c3d8fbar
654 
656  !----------------------------------------------------------------------*
657  subroutine tload_c3d8fbar &
658  (etype, nn, xx, yy, zz, tt, t0, &
659  gausses, vect, cdsys_id, coords)
660  !----------------------------------------------------------------------*
661 
662  use m_fstr
663  use mmechgauss
664  use m_matmatrix
665  use m_utilities
666 
667  !---------------------------------------------------------------------
668 
669  integer(kind=kint), parameter :: ndof = 3
670  integer(kind=kint), intent(in) :: etype, nn
671  type(tgaussstatus), intent(in) :: gausses(:)
672  real(kind=kreal), intent(in) :: xx(nn), yy(nn), zz(nn)
673  real(kind=kreal), intent(in) :: tt(nn), t0(nn)
674  real(kind=kreal), intent(out) :: vect(nn*ndof)
675  integer(kind=kint), intent(in) :: cdsys_ID
676  real(kind=kreal), intent(inout) :: coords(3, 3)
677 
678  !---------------------------------------------------------------------
679 
680  real(kind=kreal) :: alp, alp0, d(6, 6), b(6, ndof*nn)
681  real(kind=kreal) :: z1(3), det, ecoord(3, nn)
682  integer(kind=kint) :: j, LX, serr
683  real(kind=kreal) :: estrain(6), sgm(6), h(nn)
684  real(kind=kreal) :: naturalcoord(3), gderiv(nn, 3)
685  real(kind=kreal) :: wg, outa(1), ina(1), gderiv1_ave(nn, 3), alpo(3), alpo0(3), coordsys(3, 3)
686  real(kind=kreal) :: tempc, temp0, v0, jacob_ave, thermal_eps, tm(6,6)
687  logical :: ierr, matlaniso
688 
689  !---------------------------------------------------------------------
690 
691  matlaniso = .false.
692 
693  if( cdsys_id > 0 ) then ! cannot define aniso expansion when no local coord defined
694  ina = tt(1)
695  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
696  if( .not. ierr ) matlaniso = .true.
697  end if
698 
699  vect(:) = 0.0d0
700 
701  ecoord(1, :) = xx(:)
702  ecoord(2, :) = yy(:)
703  ecoord(3, :) = zz(:)
704 
705  !cal volumetric average of J=detF and dN/dx
706  v0 = 0.d0
707  jacob_ave = 0.d0
708  gderiv1_ave(1:nn,1:ndof) = 0.d0
709  do lx = 1, numofquadpoints(etype)
710  call getquadpoint( etype, lx, naturalcoord(:) )
711  call getglobalderiv( etype, nn, naturalcoord, ecoord, det, gderiv)
712  wg = getweight(etype, lx)*det
713  v0 = v0 + wg
714  jacob_ave = jacob_ave + wg
715  gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof) + wg*gderiv(1:nn, 1:ndof)
716  enddo
717  if( dabs(v0) > 1.d-12 ) then
718  if( dabs(jacob_ave) < 1.d-12 ) stop 'Error in TLOAD_C3D8Fbar: too small average jacob'
719  jacob_ave = jacob_ave/v0
720  gderiv1_ave(1:nn,1:ndof) = gderiv1_ave(1:nn,1:ndof)/(v0*jacob_ave)
721  else
722  stop 'Error in TLOAD_C3D8Fbar: too small element volume'
723  end if
724 
725  ! LOOP FOR INTEGRATION POINTS
726  do lx = 1, numofquadpoints(etype)
727 
728  call getquadpoint( etype, lx, naturalcoord(:) )
729  call getshapefunc( etype, naturalcoord, h(1:nn) )
730  call getglobalderiv(etype, nn, naturalcoord, ecoord, det, gderiv)
731 
732  if( cdsys_id > 0 ) then
733  call set_localcoordsys(coords, g_localcoordsys(cdsys_id), coordsys, serr)
734  if( serr == -1 ) stop "Fail to setup local coordinate"
735  if( serr == -2 ) then
736  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
737  end if
738  end if
739 
740  ! WEIGHT VALUE AT GAUSSIAN POINT
741  wg = getweight(etype, lx)*det
742  b(1:6, 1:nn*ndof) = 0.0d0
743  do j = 1, nn
744  b(1,3*j-2) = gderiv(j, 1)
745  b(2,3*j-1) = gderiv(j, 2)
746  b(3,3*j ) = gderiv(j, 3)
747  b(4,3*j-2) = gderiv(j, 2)
748  b(4,3*j-1) = gderiv(j, 1)
749  b(5,3*j-1) = gderiv(j, 3)
750  b(5,3*j ) = gderiv(j, 2)
751  b(6,3*j-2) = gderiv(j, 3)
752  b(6,3*j ) = gderiv(j, 1)
753  end do
754 
755  do j = 1, nn
756  z1(1:3) = (gderiv1_ave(j,1:3)-gderiv(j,1:3))/3.d0
757  b(1,3*j-2:3*j) = b(1,3*j-2:3*j)+z1(1:3)
758  b(2,3*j-2:3*j) = b(2,3*j-2:3*j)+z1(1:3)
759  b(3,3*j-2:3*j) = b(3,3*j-2:3*j)+z1(1:3)
760  end do
761 
762  tempc = dot_product( h(1:nn),tt(1:nn) )
763  temp0 = dot_product( h(1:nn),t0(1:nn) )
764 
765  call matlmatrix( gausses(lx), d3, d, 1.d0, 0.0d0, coordsys, tempc )
766 
767  ina(1) = tempc
768  if( matlaniso ) then
769  call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo(:), ierr, ina )
770  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
771  else
772  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
773  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
774  alp = outa(1)
775  end if
776  ina(1) = temp0
777  if( matlaniso ) then
778  call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo0(:), ierr, ina )
779  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
780  else
781  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
782  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
783  alp0 = outa(1)
784  end if
785 
786  !**
787  !** THERMAL strain
788  !**
789  if( matlaniso ) then
790  do j = 1, 3
791  estrain(j) = alpo(j)*(tempc-ref_temp)-alpo0(j)*(temp0-ref_temp)
792  end do
793  estrain(4:6) = 0.0d0
794  call transformation(coordsys, tm)
795  estrain(:) = matmul( estrain(:), tm ) ! to global coord
796  estrain(4:6) = estrain(4:6)*2.0d0
797  else
798  thermal_eps = alp*(tempc-ref_temp)-alp0*(temp0-ref_temp)
799  estrain(1:3) = thermal_eps
800  estrain(4:6) = 0.0d0
801  end if
802 
803  !**
804  !** SET SGM {s}=[D]{e}
805  !**
806  sgm(:) = matmul( d(:, :), estrain(:) )
807 
808  !**
809  !** CALCULATE LOAD {F}=[B]T{e}
810  !**
811  vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(:), b(:, :) )*wg
812 
813  end do
814 
815  end subroutine tload_c3d8fbar
816 
817 
818 end module m_static_lib_fbar
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
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_fbar::tload_c3d8fbar
subroutine tload_c3d8fbar(etype, nn, XX, YY, ZZ, TT, T0, gausses, VECT, cdsys_ID, coords)
This subroutien calculate thermal loading.
Definition: static_LIB_Fbar.f90:660
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_static_lib_fbar::stf_c3d8fbar
subroutine stf_c3d8fbar(etype, nn, ecoord, gausses, stiff, cdsys_ID, coords, time, tincr, u, temperature)
This subroutine calculate stiff matrix using F-bar method.
Definition: static_LIB_Fbar.f90:25
m_common_struct::g_localcoordsys
type(tlocalcoordsys), dimension(:), pointer, save g_localcoordsys
Definition: m_common_struct.f90:19
m_static_lib_fbar
This module contains several strategy to free locking problem in Eight-node hexagonal element.
Definition: static_LIB_Fbar.f90:7
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_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
m_static_lib_fbar::update_c3d8fbar
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.
Definition: static_LIB_Fbar.f90:340
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_utilities::determinant33
real(kind=kreal) function determinant33(XJ)
Compute determinant for 3*3 matrix.
Definition: utilities.f90:233
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