FrontISTR  5.7.1
Large-scale structural analysis program with finit element method
static_LIB_3dIC.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 !-------------------------------------------------------------------------------
9 module m_static_lib_3dic
10 
11  use hecmw, only : kint, kreal
12  use m_utilities
13  use elementinfo
14 
15  implicit none
16 
17 contains
18 
20  !----------------------------------------------------------------------*
21  subroutine stf_c3d8ic &
22  (etype, nn, ecoord, gausses, stiff, cdsys_id, coords, &
23  time, tincr, u, aux, temperature)
24  !----------------------------------------------------------------------*
25 
26  use mmechgauss
27  use m_matmatrix
28  use m_common_struct
29  use m_static_lib_3d, only: geomat_c3
30 
31  !---------------------------------------------------------------------
32 
33  integer(kind=kint), intent(in) :: etype
34  integer(kind=kint), intent(in) :: nn
35  real(kind=kreal), intent(in) :: ecoord(3, nn)
36  type(tgaussstatus), intent(in) :: gausses(:)
37  real(kind=kreal), intent(out) :: stiff(:, :)
38  integer(kind=kint), intent(in) :: cdsys_id
39  real(kind=kreal), intent(inout) :: coords(3, 3)
40  real(kind=kreal), intent(in) :: time
41  real(kind=kreal), intent(in) :: tincr
42  real(kind=kreal), intent(in), optional :: u(3, nn)
43  real(kind=kreal), intent(in), optional :: aux(3, 3)
44  real(kind=kreal), intent(in) :: temperature(nn)
45 
46  !---------------------------------------------------------------------
47 
48  integer(kind=kint) :: flag
49  integer(kind=kint), parameter :: ndof = 3
50  real(kind=kreal) :: d(6, 6), b(6, ndof*(nn+3)), db(6, ndof*(nn+3))
51  real(kind=kreal) :: gderiv(nn+3, 3), stress(6)
52  real(kind=kreal) :: xj(9, 9), jacobian(3, 3), inverse(3, 3)
53  real(kind=kreal) :: tmpstiff((nn+3)*3, (nn+3)*3), tmpk(nn*3, 9)
54  real(kind=kreal) :: det, wg, elem(3, nn), mat(6, 6)
55  integer(kind=kint) :: i, j, lx
56  integer(kind=kint) :: serr
57  real(kind=kreal) :: naturalcoord(3), unode(3, nn+3)
58  real(kind=kreal) :: gdispderiv(3, 3), coordsys(3, 3)
59  real(kind=kreal) :: b1(6, ndof*(nn+3))
60  real(kind=kreal) :: smat(9, 9)
61  real(kind=kreal) :: bn(9, ndof*(nn+3)), sbn(9, ndof*(nn+3))
62  real(kind=kreal) :: spfunc(nn), temp
63 
64  !---------------------------------------------------------------------
65 
66  if( present(u) .AND. present(aux) ) then
67  unode(:, 1:nn) = u(:, :)
68  unode(:, nn+1:nn+3) = aux(:, :)
69  end if
70 
71  ! we suppose the same material type in the element
72  flag = gausses(1)%pMaterial%nlgeom_flag
73  if( .not. present(u) ) flag = infinitesimal ! enforce to infinitesimal deformation analysis
74  elem(:, :) = ecoord(:, :)
75  if( flag == updatelag ) elem(:, :) = ecoord(:, :)+unode(:, 1:nn)
76 
77  ! --- Inverse of Jacobian at elemental center
78  naturalcoord(:) = 0.0d0
79  call getjacobian(etype, nn, naturalcoord, elem, det, jacobian, inverse)
80  inverse(:, :)= inverse(:, :)*det
81  ! ---- We now calculate stiff matrix include incompatible mode
82  ! [ Kdd Kda ]
83  ! [ Kad Kaa ]
84  tmpstiff(:, :) = 0.0d0
85  b1(1:6, 1:(nn+3)*ndof) = 0.0d0
86  bn(1:9, 1:(nn+3)*ndof) = 0.0d0
87 
88  do lx = 1, numofquadpoints(etype)
89 
90  call getquadpoint(etype, lx, naturalcoord)
91  call getglobalderiv( etype, nn, naturalcoord, elem, det, gderiv(1:nn, 1:3) )
92 
93  if( cdsys_id > 0 ) then
94  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
95  if( serr == -1 ) stop "Fail to setup local coordinate"
96  if( serr == -2 ) then
97  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
98  end if
99  end if
100 
101  call getshapefunc( etype, naturalcoord, spfunc )
102  temp = dot_product( temperature, spfunc )
103  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
104 
105  if( flag == updatelag ) then
106  call geomat_c3( gausses(lx)%stress, mat )
107  d(:, :) = d(:, :)-mat
108  endif
109 
110  ! -- Derivative of shape function of incompatible mode --
111  ! [ -2*a 0, 0 ]
112  ! [ 0, -2*b, 0 ]
113  ! [ 0, 0, -2*c ]
114  ! we don't call getShapeDeriv but use above value directly for computation efficiency
115  gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse(1, :)/det
116  gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse(2, :)/det
117  gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse(3, :)/det
118 
119  wg = getweight(etype, lx)*det
120  b(1:6, 1:(nn+3)*ndof)=0.0d0
121  do j = 1, nn+3
122  b(1, 3*j-2) = gderiv(j, 1)
123  b(2, 3*j-1) = gderiv(j, 2)
124  b(3, 3*j ) = gderiv(j, 3)
125  b(4, 3*j-2) = gderiv(j, 2)
126  b(4, 3*j-1) = gderiv(j, 1)
127  b(5, 3*j-1) = gderiv(j, 3)
128  b(5, 3*j ) = gderiv(j, 2)
129  b(6, 3*j-2) = gderiv(j, 3)
130  b(6, 3*j ) = gderiv(j, 1)
131  end do
132 
133  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
134  if( flag == totallag ) then
135  ! ---dudx(i,j) ==> gdispderiv(i,j)
136  gdispderiv(1:ndof,1:ndof) = matmul( unode(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) )
137  do j = 1, nn+3
138 
139  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
140  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
141  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
142  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
143  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
144  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
145  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
146  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
147  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
148  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
149  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
150  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
151  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
152  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
153  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
154  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
155  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
156  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
157 
158  end do
159  ! ---BL = BL0 + BL1
160  do j = 1, (nn+3)*ndof
161  b(:, j) = b(:, j)+b1(:, j)
162  end do
163 
164  end if
165 
166  db(1:6, 1:(nn+3)*ndof) = matmul( d, b(1:6, 1:(nn+3)*ndof) )
167  do j=1,(nn+3)*ndof
168  do i=1,(nn+3)*ndof
169  tmpstiff(i, j) = tmpstiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
170  enddo
171  enddo
172 
173  ! calculate the stress matrix ( TOTAL LAGRANGE METHOD )
174  if( flag == totallag .OR. flag == updatelag ) then
175  stress(1:6) = gausses(lx)%stress
176  do j = 1, nn+3
177  bn(1, 3*j-2) = gderiv(j, 1)
178  bn(2, 3*j-1) = gderiv(j, 1)
179  bn(3, 3*j ) = gderiv(j, 1)
180  bn(4, 3*j-2) = gderiv(j, 2)
181  bn(5, 3*j-1) = gderiv(j, 2)
182  bn(6, 3*j ) = gderiv(j, 2)
183  bn(7, 3*j-2) = gderiv(j, 3)
184  bn(8, 3*j-1) = gderiv(j, 3)
185  bn(9, 3*j ) = gderiv(j, 3)
186  end do
187  smat(:, :) = 0.0d0
188  do j = 1, 3
189  smat(j , j ) = stress(1)
190  smat(j , j+3) = stress(4)
191  smat(j , j+6) = stress(6)
192  smat(j+3, j ) = stress(4)
193  smat(j+3, j+3) = stress(2)
194  smat(j+3, j+6) = stress(5)
195  smat(j+6, j ) = stress(6)
196  smat(j+6, j+3) = stress(5)
197  smat(j+6, j+6) = stress(3)
198  end do
199  sbn(1:9, 1:(nn+3)*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:(nn+3)*ndof) )
200  do j=1,(nn+3)*ndof
201  do i=1,(nn+3)*ndof
202  tmpstiff(i, j) = tmpstiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
203  enddo
204  enddo
205  end if
206 
207  end do
208 
209  ! -----Condense tmpstiff to stiff
210  xj(1:9, 1:9)= tmpstiff(nn*ndof+1:(nn+3)*ndof, nn*ndof+1:(nn+3)*ndof)
211  call calinverse(9, xj)
212  tmpk = matmul( tmpstiff( 1:nn*ndof, nn*ndof+1:(nn+3)*ndof ), xj )
213  stiff(1:nn*ndof, 1:nn*ndof) = tmpstiff(1:nn*ndof, 1:nn*ndof)-matmul( tmpk, tmpstiff(nn*ndof+1:(nn+3)*ndof, 1:nn*ndof) )
214 
215  end subroutine stf_c3d8ic
216 
217 
219  !---------------------------------------------------------------------*
220  subroutine update_c3d8ic &
221  (etype,nn,ecoord, u, du, ddu, cdsys_id, coords, &
222  qf, gausses, iter, time, tincr, aux, ddaux, tt, t0, tn )
223  !---------------------------------------------------------------------*
224 
225  use m_fstr
226  use mmaterial
227  use mmechgauss
228  use m_matmatrix
229  use m_elastoplastic
230  use m_utilities
231  use m_static_lib_3d
232 
233  integer(kind=kint), intent(in) :: etype
234  integer(kind=kint), intent(in) :: nn
235  real(kind=kreal), intent(in) :: ecoord(3, nn)
236  real(kind=kreal), intent(in) :: u(3, nn)
237  real(kind=kreal), intent(in) :: du(3, nn)
238  real(kind=kreal), intent(in) :: ddu(3, nn)
239  integer(kind=kint), intent(in) :: cdsys_id
240  real(kind=kreal), intent(inout) :: coords(3, 3)
241  real(kind=kreal), intent(out) :: qf(nn*3)
242  type(tgaussstatus), intent(inout) :: gausses(:)
243  integer, intent(in) :: iter
244  real(kind=kreal), intent(in) :: time
245  real(kind=kreal), intent(in) :: tincr
246  real(kind=kreal), intent(inout) :: aux(3, 3)
247  real(kind=kreal), intent(out) :: ddaux(3, 3)
248  real(kind=kreal), intent(in) :: tt(nn)
249  real(kind=kreal), intent(in) :: t0(nn)
250  real(kind=kreal), intent(in) :: tn(nn)
251 
252  ! LOCAL VARIABLES
253  integer(kind=kint) :: flag
254  integer(kind=kint), parameter :: ndof = 3
255  real(kind=kreal) :: d(6,6), b(6,ndof*(nn+3)), b1(6,ndof*(nn+3)), spfunc(nn), ina(1)
256  real(kind=kreal) :: bn(9,ndof*(nn+3)), db(6, ndof*(nn+3)), stress(6), smat(9, 9), sbn(9, ndof*(nn+3))
257  real(kind=kreal) :: gderiv(nn+3,3), gderiv0(nn+3,3), gdispderiv(3,3), f(3,3), det, det0, wg, ttc, tt0, ttn
258  integer(kind=kint) :: i, j, lx, mtype, serr
259  real(kind=kreal) :: naturalcoord(3), rot(3,3), mat(6,6), epsth(6)
260  real(kind=kreal) :: totaldisp(3,nn+3), elem(3,nn), elem1(3,nn), coordsys(3,3)
261  real(kind=kreal) :: dstrain(6)
262  real(kind=kreal) :: alpo(3)
263  logical :: ierr, matlaniso
264  real(kind=kreal) :: stiff_ad(9, nn*3), stiff_aa(9, 9), qf_a(9)
265  real(kind=kreal) :: xj(9, 9)
266  real(kind=kreal) :: tmpforce(9), tmpdisp(9), tmp_d(nn*3), tmp_a(9)
267  real(kind=kreal) :: jacobian(3, 3), inverse(3, 3), inverse1(3, 3), inverse0(3, 3)
268 
269  !---------------------------------------------------------------------
270 
271  totaldisp(:, 1:nn) = u(:, :) + du(:, :) - ddu(:, :)
272  totaldisp(:, nn+1:nn+3) = aux(:, :)
273 
274  ! we suppose the same material type in the element
275  flag = gausses(1)%pMaterial%nlgeom_flag
276  elem(:, :) = ecoord(:, :)
277  if( flag == updatelag ) elem(:, :) = ecoord(:, :)+totaldisp(:, 1:nn)
278 
279  ! --- Inverse of Jacobian at elemental center
280  naturalcoord(:) = 0.0d0
281  call getjacobian(etype, nn, naturalcoord, elem, det, jacobian, inverse)
282  inverse(:, :)= inverse(:, :)*det
283  ! ---- We now calculate stiff matrix include incompatible mode
284  ! [ Kdd Kda ]
285  ! [ Kad Kaa ]
286  stiff_ad(:, :) = 0.0d0
287  stiff_aa(:, :) = 0.0d0
288  b1(1:6, 1:(nn+3)*ndof) = 0.0d0
289  bn(1:9, 1:(nn+3)*ndof) = 0.0d0
290 
291  do lx = 1, numofquadpoints(etype)
292 
293  call getquadpoint(etype, lx, naturalcoord)
294  call getglobalderiv( etype, nn, naturalcoord, elem, det, gderiv(1:nn, 1:3) )
295 
296  if( cdsys_id > 0 ) then
297  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
298  if( serr == -1 ) stop "Fail to setup local coordinate"
299  if( serr == -2 ) then
300  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
301  end if
302  end if
303 
304  call getshapefunc( etype, naturalcoord, spfunc )
305  ttc = dot_product( tt, spfunc )
306  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, ttc )
307 
308  if( flag == updatelag ) then
309  call geomat_c3( gausses(lx)%stress, mat )
310  d(:, :) = d(:, :)-mat
311  endif
312 
313  ! -- Derivative of shape function of incompatible mode --
314  ! [ -2*a 0, 0 ]
315  ! [ 0, -2*b, 0 ]
316  ! [ 0, 0, -2*c ]
317  ! we don't call getShapeDeriv but use above value directly for computation efficiency
318  gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse(1, :)/det
319  gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse(2, :)/det
320  gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse(3, :)/det
321 
322  wg = getweight(etype, lx)*det
323  b(1:6, 1:(nn+3)*ndof)=0.0d0
324  do j = 1, nn+3
325  b(1, 3*j-2) = gderiv(j, 1)
326  b(2, 3*j-1) = gderiv(j, 2)
327  b(3, 3*j ) = gderiv(j, 3)
328  b(4, 3*j-2) = gderiv(j, 2)
329  b(4, 3*j-1) = gderiv(j, 1)
330  b(5, 3*j-1) = gderiv(j, 3)
331  b(5, 3*j ) = gderiv(j, 2)
332  b(6, 3*j-2) = gderiv(j, 3)
333  b(6, 3*j ) = gderiv(j, 1)
334  end do
335 
336  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
337  if( flag == totallag ) then
338  ! ---dudx(i,j) ==> gdispderiv(i,j)
339  gdispderiv(1:ndof,1:ndof) = matmul( totaldisp(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) )
340  do j = 1, nn+3
341 
342  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
343  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
344  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
345  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
346  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
347  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
348  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
349  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
350  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
351  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
352  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
353  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
354  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
355  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
356  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
357  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
358  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
359  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
360 
361  end do
362  ! ---BL = BL0 + BL1
363  do j = 1, (nn+3)*ndof
364  b(:, j) = b(:, j)+b1(:, j)
365  end do
366 
367  end if
368 
369  db(1:6, 1:(nn+3)*ndof) = matmul( d, b(1:6, 1:(nn+3)*ndof) )
370  do j=1,nn*ndof
371  do i=1,3*ndof
372  stiff_ad(i, j) = stiff_ad(i, j)+dot_product( b(:, nn*ndof+i), db(:, j) )*wg
373  enddo
374  enddo
375 
376  do j=1,3*ndof
377  do i=1,3*ndof
378  stiff_aa(i, j) = stiff_aa(i, j)+dot_product( b(:, nn*ndof+i), db(:, nn*ndof+j) )*wg
379  enddo
380  enddo
381 
382  ! calculate the stress matrix ( TOTAL LAGRANGE METHOD )
383  if( flag == totallag .OR. flag == updatelag ) then
384  stress(1:6) = gausses(lx)%stress
385  do j = 1, nn+3
386  bn(1, 3*j-2) = gderiv(j, 1)
387  bn(2, 3*j-1) = gderiv(j, 1)
388  bn(3, 3*j ) = gderiv(j, 1)
389  bn(4, 3*j-2) = gderiv(j, 2)
390  bn(5, 3*j-1) = gderiv(j, 2)
391  bn(6, 3*j ) = gderiv(j, 2)
392  bn(7, 3*j-2) = gderiv(j, 3)
393  bn(8, 3*j-1) = gderiv(j, 3)
394  bn(9, 3*j ) = gderiv(j, 3)
395  end do
396  smat(:, :) = 0.0d0
397  do j = 1, 3
398  smat(j , j ) = stress(1)
399  smat(j , j+3) = stress(4)
400  smat(j , j+6) = stress(6)
401  smat(j+3, j ) = stress(4)
402  smat(j+3, j+3) = stress(2)
403  smat(j+3, j+6) = stress(5)
404  smat(j+6, j ) = stress(6)
405  smat(j+6, j+3) = stress(5)
406  smat(j+6, j+6) = stress(3)
407  end do
408  sbn(1:9, 1:(nn+3)*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:(nn+3)*ndof) )
409 
410  do j=1,nn*ndof
411  do i=1,3*ndof
412  stiff_ad(i, j) = stiff_ad(i, j)+dot_product( bn(:, nn*ndof+i), sbn(:, j) )*wg
413  enddo
414  enddo
415  do j=1,3*ndof
416  do i=1,3*ndof
417  stiff_aa(i, j) = stiff_aa(i, j)+dot_product( bn(:, nn*ndof+i), sbn(:, nn*ndof+j) )*wg
418  enddo
419  enddo
420  end if
421 
422  end do
423 
424  ! -----recover incompatible dof of nodal displacement
425  xj(1:3*ndof, 1:3*ndof)= stiff_aa(1:3*ndof, 1:3*ndof)
426  call calinverse(3*ndof, xj)
427  ! --- [Kad] * ddu
428  do i=1,nn
429  tmp_d((i-1)*ndof+1:i*ndof) = ddu(1:ndof, i)
430  enddo
431  tmpforce(1:3*ndof) = matmul( stiff_ad(1:3*ndof, 1:nn*ndof), tmp_d(1:nn*ndof) )
432  ! --- ddaux = -[Kaa]-1 * ([Kad] * ddu)
433  tmp_a(1:3*ndof) = -matmul( xj(1:3*ndof, 1:3*ndof), tmpforce(1:3*ndof) )
434  do i=1,3
435  ddaux(1:ndof, i) = tmp_a((i-1)*ndof+1:i*ndof)
436  enddo
437 
438 
439  !---------------------------------------------------------------------
440 
441  totaldisp(:,1:nn) = u(:,:)+ du(:,:)
442  totaldisp(:,nn+1:nn+3) = aux(:,:) + ddaux(:,:)
443 
444  !---------------------------------------------------------------------
445 
446  qf(:) = 0.0d0
447  qf_a(:) = 0.0d0
448 
449  stiff_ad(:, :) = 0.0d0
450  stiff_aa(:, :) = 0.0d0
451 
452  !---------------------------------------------------------------------
453 
454  if( flag == updatelag ) then
455  elem(:,:) = (0.5d0*du(:,:)+u(:,:) ) +ecoord(:,:)
456  elem1(:,:) = (du(:,:)+u(:,:) ) +ecoord(:,:)
457  ! elem = elem1
458  totaldisp(:,1:nn) = du(:,:)
459  if( iter == 1 ) aux(:,:) = 0.0d0 !--- is this correct???
460  totaldisp(:,nn+1:nn+3) = aux(:,:)
461  end if
462 
463  matlaniso = .false.
464  ina = tt(1)
465  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
466  if( .not. ierr ) matlaniso = .true.
467 
468  ! --- Inverse of Jacobian at elemental center
469  naturalcoord(:) = 0.0d0
470  call getjacobian(etype, nn, naturalcoord, elem, det, jacobian, inverse)
471  inverse(:, :)= inverse(:, :)*det
472  if( flag == updatelag ) then
473  call getjacobian(etype, nn, naturalcoord, elem1, det, jacobian, inverse1)
474  inverse1(:, :)= inverse1(:, :)*det
475  call getjacobian(etype, nn, naturalcoord, ecoord, det, jacobian, inverse0)
476  inverse0(:, :)= inverse0(:, :)*det !for deformation gradient F
477  endif
478 
479  do lx = 1, numofquadpoints(etype)
480 
481  mtype = gausses(lx)%pMaterial%mtype
482 
483  call getquadpoint( etype, lx, naturalcoord(:) )
484  call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
485 
486  if( cdsys_id > 0 ) then
487  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:,:), serr )
488  if( serr == -1 ) stop "Fail to setup local coordinate"
489  if( serr == -2 ) then
490  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
491  end if
492  end if
493 
494  ! ========================================================
495  ! UPDATE STRAIN and STRESS
496  ! ========================================================
497 
498  ! Thermal Strain
499  epsth = 0.0d0
500  call getshapefunc(etype, naturalcoord, spfunc)
501  ttc = dot_product(tt, spfunc)
502  tt0 = dot_product(t0, spfunc)
503  ttn = dot_product(tn, spfunc)
504  call cal_thermal_expansion_c3( tt0, ttc, gausses(lx)%pMaterial, coordsys, matlaniso, epsth )
505 
506  ! -- Derivative of shape function of incompatible mode --
507  ! [ -2*a 0, 0 ]
508  ! [ 0, -2*b, 0 ]
509  ! [ 0, 0, -2*c ]
510  ! we don't call getShapeDeriv but use above value directly for computation efficiency
511  gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse(1, :)/det
512  gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse(2, :)/det
513  gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse(3, :)/det
514 
515  ! Small strain
516  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) )
517  dstrain(1) = gdispderiv(1, 1)
518  dstrain(2) = gdispderiv(2, 2)
519  dstrain(3) = gdispderiv(3, 3)
520  dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
521  dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
522  dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
523  dstrain(:) = dstrain(:)-epsth(:) ! alright?
524 
525  f(1:3,1:3) = 0.d0; f(1,1)=1.d0; f(2,2)=1.d0; f(3,3)=1.d0; !deformation gradient
526  if( flag == infinitesimal ) then
527  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
528  f(1:3,1:3) = f(1:3,1:3) + gdispderiv(1:3,1:3)
529 
530  else if( flag == totallag ) then
531  ! Green-Lagrange strain
532  dstrain(1) = dstrain(1)+0.5d0*dot_product( gdispderiv(:, 1), gdispderiv(:, 1) )
533  dstrain(2) = dstrain(2)+0.5d0*dot_product( gdispderiv(:, 2), gdispderiv(:, 2) )
534  dstrain(3) = dstrain(3)+0.5d0*dot_product( gdispderiv(:, 3), gdispderiv(:, 3) )
535  dstrain(4) = dstrain(4)+( gdispderiv(1, 1)*gdispderiv(1, 2) &
536  +gdispderiv(2, 1)*gdispderiv(2, 2)+gdispderiv(3, 1)*gdispderiv(3, 2) )
537  dstrain(5) = dstrain(5)+( gdispderiv(1, 2)*gdispderiv(1, 3) &
538  +gdispderiv(2, 2)*gdispderiv(2, 3)+gdispderiv(3, 2)*gdispderiv(3, 3) )
539  dstrain(6) = dstrain(6)+( gdispderiv(1, 1)*gdispderiv(1, 3) &
540  +gdispderiv(2, 1)*gdispderiv(2, 3)+gdispderiv(3, 1)*gdispderiv(3, 3) )
541 
542  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
543 
544  else if( flag == updatelag ) then
545  ! CALL GEOMAT_C3( gausses(LX)%stress, mat )
546  ! D(:, :) = D(:, :)+mat(:, :)
547  rot = 0.0d0
548  rot(1, 2)= 0.5d0*(gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
549  rot(2, 3)= 0.5d0*(gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
550  rot(1, 3)= 0.5d0*(gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
551 
552  gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+ dstrain(1:6)+epsth(:)
553  call getglobalderiv(etype, nn, naturalcoord, ecoord, det0, gderiv0)
554  gderiv0(nn+1, :) = -2.0d0*naturalcoord(1)*inverse0(1, :)/det0
555  gderiv0(nn+2, :) = -2.0d0*naturalcoord(2)*inverse0(2, :)/det0
556  gderiv0(nn+3, :) = -2.0d0*naturalcoord(3)*inverse0(3, :)/det0
557  f(1:3,1:3) = f(1:3,1:3) + matmul( totaldisp(1:ndof, 1:nn+3), gderiv0(1:nn+3, 1:ndof) )
558 
559  end if
560 
561  ! Update stress
562  call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn )
563 
564  ! ========================================================
565  ! calculate the internal force ( equivalent nodal force )
566  ! ========================================================
567  ! Small strain
568  b(1:6, 1:(nn+3)*ndof) = 0.0d0
569  do j=1,nn+3
570  b(1,3*j-2) = gderiv(j, 1)
571  b(2,3*j-1) = gderiv(j, 2)
572  b(3,3*j ) = gderiv(j, 3)
573  b(4,3*j-2) = gderiv(j, 2)
574  b(4,3*j-1) = gderiv(j, 1)
575  b(5,3*j-1) = gderiv(j, 3)
576  b(5,3*j ) = gderiv(j, 2)
577  b(6,3*j-2) = gderiv(j, 3)
578  b(6,3*j ) = gderiv(j, 1)
579  end do
580 
581  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
582  if( flag == infinitesimal ) then
583 
584  else if( flag == totallag ) then
585 
586  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) )
587  b1(1:6, 1:(nn+3)*ndof)=0.0d0
588  do j = 1,nn+3
589  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
590  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
591  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
592  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
593  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
594  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
595  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
596  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
597  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
598  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
599  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
600  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
601  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
602  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
603  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
604  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
605  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
606  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
607  end do
608  ! BL = BL0 + BL1
609  do j=1,(nn+3)*ndof
610  b(:,j) = b(:,j)+b1(:,j)
611  end do
612 
613  else if( flag == updatelag ) then
614 
615  call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv(1:nn,1:3))
616 
617  ! -- Derivative of shape function of incompatible mode --
618  ! [ -2*a 0, 0 ]
619  ! [ 0, -2*b, 0 ]
620  ! [ 0, 0, -2*c ]
621  ! we don't call getShapeDeriv but use above value directly for computation efficiency
622  gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse1(1, :)/det
623  gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse1(2, :)/det
624  gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse1(3, :)/det
625 
626  b(1:6, 1:(nn+3)*ndof) = 0.0d0
627  do j = 1, nn+3
628  b(1, 3*j-2) = gderiv(j, 1)
629  b(2, 3*j-1) = gderiv(j, 2)
630  b(3, 3*j ) = gderiv(j, 3)
631  b(4, 3*j-2) = gderiv(j, 2)
632  b(4, 3*j-1) = gderiv(j, 1)
633  b(5, 3*j-1) = gderiv(j, 3)
634  b(5, 3*j ) = gderiv(j, 2)
635  b(6, 3*j-2) = gderiv(j, 3)
636  b(6, 3*j ) = gderiv(j, 1)
637  end do
638 
639  end if
640 
641  wg=getweight( etype, lx )*det
642 
643  db(1:6, 1:(nn+3)*ndof) = matmul( d, b(1:6, 1:(nn+3)*ndof) )
644 
645  do j=1,nn*ndof
646  do i=1,3*ndof
647  stiff_ad(i, j) = stiff_ad(i, j)+dot_product( b(:, nn*ndof+i), db(:, j) )*wg
648  enddo
649  enddo
650  do j=1,3*ndof
651  do i=1,3*ndof
652  stiff_aa(i, j) = stiff_aa(i, j)+dot_product( b(:, nn*ndof+i), db(:, nn*ndof+j) )*wg
653  enddo
654  enddo
655 
656  ! calculate the Internal Force
657  qf(1:nn*ndof) &
658  = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
659  qf_a(1:3*ndof) &
660  = qf_a(1:3*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,nn*ndof+1:(nn+3)*ndof) )*wg
661 
662  ! integrate strain energy
663  gausses(lx)%strain_energy = gausses(lx)%strain_energy*wg
664  end do
665 
666  ! condence ( qf - [Kda] * [Kaa]-1 * qf_a )
667  xj(1:9, 1:9)= stiff_aa(1:3*ndof, 1:3*ndof)
668  call calinverse(9, xj)
669  tmpdisp(:) = matmul( xj(:,:), qf_a(:) )
670  do i=1,nn*ndof
671  qf(i) = qf(i) - dot_product( stiff_ad(:,i), tmpdisp(:) )
672  enddo
673  end subroutine update_c3d8ic
674 
675 
677  !----------------------------------------------------------------------*
678  subroutine tload_c3d8ic &
679  (etype, nn, xx, yy, zz, tt, t0, &
680  gausses, vect, cdsys_id, coords)
681  !----------------------------------------------------------------------*
682 
683  use m_fstr
684  use mmechgauss
685  use m_matmatrix
686  use m_common_struct
687 
688  !---------------------------------------------------------------------
689 
690  integer(kind=kint), parameter :: ndof=3
691  integer(kind=kint), intent(in) :: etype
692  integer(kind=kint), intent(in) :: nn
693  real(kind=kreal), intent(in) :: xx(nn), yy(nn), zz(nn)
694  real(kind=kreal), intent(in) :: tt(nn), t0(nn)
695  type(tgaussstatus), intent(inout) :: gausses(:)
696  real(kind=kreal), intent(out) :: vect(nn*ndof)
697  integer(kind=kint), intent(in) :: cdsys_id
698  real(kind=kreal), intent(inout) :: coords(3, 3)
699 
700  !---------------------------------------------------------------------
701 
702  real(kind=kreal) :: alp, alp0
703  real(kind=kreal) :: d(6, 6), b(6, ndof*(nn+3)), db_a(6, ndof*3)
704  real(kind=kreal) :: det, wg, ecoord(3, nn)
705  integer(kind=kint) :: i, j, ic, serr
706  real(kind=kreal) :: epsth(6), sgm(6), h(nn), alpo(3), alpo0(3), coordsys(3, 3), xj(9, 9)
707  real(kind=kreal) :: naturalcoord(3), gderiv(nn+3, 3)
708  real(kind=kreal) :: jacobian(3, 3),inverse(3, 3)
709  real(kind=kreal) :: stiff_da(nn*3, 3*3), stiff_aa(3*3, 3*3)
710  real(kind=kreal) :: vect_a(3*3)
711  real(kind=kreal) :: icdisp(9)
712  real(kind=kreal) :: tempc, temp0, thermal_eps, tm(6, 6), outa(1), ina(1)
713  logical :: ierr, matlaniso
714 
715  !---------------------------------------------------------------------
716 
717  matlaniso = .false.
718  if( cdsys_id > 0 ) then ! cannot define aniso expansion when no local coord defined
719  ina = tt(1)
720  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
721  if( .not. ierr ) matlaniso = .true.
722  end if
723 
724  ecoord(1, :) = xx(:)
725  ecoord(2, :) = yy(:)
726  ecoord(3, :) = zz(:)
727  ! ---- Calculate enhanced displacement at first
728  ! --- Inverse of Jacobian at elemental center
729  naturalcoord(:) = 0.0d0
730  call getjacobian(etype, nn, naturalcoord, ecoord, det, jacobian, inverse)
731  inverse(:, :) = inverse(:, :)*det
732  ! ---- We now calculate stiff matrix include incompatible mode
733  ! [ Kdd Kda ]
734  ! [ Kad Kaa ]
735  stiff_da(:, :) = 0.0d0
736  stiff_aa(:, :) = 0.0d0
737  b(1:6, 1:(nn+3)*ndof) = 0.0d0
738  vect(1:nn*ndof) = 0.0d0
739  vect_a(1:3*ndof) = 0.0d0
740 
741  do ic = 1, numofquadpoints(etype)
742 
743  call getquadpoint(etype, ic, naturalcoord)
744  call getglobalderiv( etype, nn, naturalcoord, ecoord, det, gderiv(1:nn, 1:3) )
745 
746  if( cdsys_id > 0 ) then
747  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
748  if( serr == -1 ) stop "Fail to setup local coordinate"
749  if( serr == -2 ) then
750  write(*,*) "WARNING! Cannot setup local coordinate, it is modified automatically"
751  end if
752  end if
753 
754  call getshapefunc( etype, naturalcoord, h(1:nn) )
755  tempc = dot_product( h(1:nn), tt(1:nn) )
756  temp0 = dot_product( h(1:nn), t0(1:nn) )
757  call matlmatrix( gausses(ic), d3, d, 1.d0, 1.0d0, coordsys, tempc )
758 
759  ! -- Derivative of shape function of incompatible mode --
760  ! [ -2*a 0, 0 ]
761  ! [ 0, -2*b, 0 ]
762  ! [ 0, 0, -2*c ]
763  ! we don't call getShapeDeriv but use above value directly for computation efficiency
764  gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse(1, :)/det
765  gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse(2, :)/det
766  gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse(3, :)/det
767 
768  wg = getweight(etype, ic)*det
769  do j = 1, nn+3
770  b(1, 3*j-2) = gderiv(j, 1)
771  b(2, 3*j-1) = gderiv(j, 2)
772  b(3, 3*j ) = gderiv(j, 3)
773  b(4, 3*j-2) = gderiv(j, 2)
774  b(4, 3*j-1) = gderiv(j, 1)
775  b(5, 3*j-1) = gderiv(j, 3)
776  b(5, 3*j ) = gderiv(j, 2)
777  b(6, 3*j-2) = gderiv(j, 3)
778  b(6, 3*j ) = gderiv(j, 1)
779  end do
780 
781  db_a(1:6, 1:3*ndof) = matmul( d, b(1:6, nn*ndof+1:(nn+3)*ndof) )
782  do j=1,3*ndof
783  do i=1,nn*ndof
784  stiff_da(i, j) = stiff_da(i, j) + dot_product( b(1:6, i), db_a(1:6, j) ) * wg
785  enddo
786  enddo
787  do j=1,3*ndof
788  do i=1,3*ndof
789  stiff_aa(i, j) = stiff_aa(i, j) + dot_product( b(1:6, nn*ndof+i), db_a(1:6, j) ) * wg
790  enddo
791  enddo
792 
793  ina(1) = tempc
794  if( matlaniso ) then
795  call fetch_tabledata( mc_orthoexp, gausses(ic)%pMaterial%dict, alpo(:), ierr, ina )
796  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
797  else
798  call fetch_tabledata( mc_themoexp, gausses(ic)%pMaterial%dict, outa(:), ierr, ina )
799  if( ierr ) outa(1) = gausses(ic)%pMaterial%variables(m_exapnsion)
800  alp = outa(1)
801  end if
802  ina(1) = temp0
803  if( matlaniso ) then
804  call fetch_tabledata( mc_orthoexp, gausses(ic)%pMaterial%dict, alpo0(:), ierr, ina )
805  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
806  else
807  call fetch_tabledata( mc_themoexp, gausses(ic)%pMaterial%dict, outa(:), ierr, ina )
808  if( ierr ) outa(1) = gausses(ic)%pMaterial%variables(m_exapnsion)
809  alp0 = outa(1)
810  end if
811 
812  !**
813  !** THERMAL strain
814  !**
815  if( matlaniso ) then
816  do j = 1, 3
817  epsth(j) = alpo(j)*(tempc-ref_temp)-alpo0(j)*(temp0-ref_temp)
818  end do
819  epsth(4:6) = 0.0d0
820  call transformation(coordsys, tm)
821  epsth(:) = matmul( epsth(:), tm ) ! to global coord
822  epsth(4:6) = epsth(4:6)*2.0d0
823  else
824  thermal_eps = alp*(tempc-ref_temp)-alp0*(temp0-ref_temp)
825  epsth(1:3) = thermal_eps
826  epsth(4:6) = 0.0d0
827  end if
828 
829  !**
830  !** SET SGM {s}=[D]{e}
831  !**
832  sgm(:) = matmul( d(:, :), epsth(:) )
833 
834  !**
835  !** CALCULATE LOAD {F}=[B]T{e}
836  !**
837  vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(1:6), b(1:6, 1:nn*ndof) )*wg
838  vect_a(1:3*ndof) = vect_a(1:3*ndof)+matmul( sgm(1:6), b(1:6, nn*ndof+1:(nn+3)*ndof) )*wg
839 
840  end do
841 
842  ! --- condense vect
843  xj(1:9,1:9)= stiff_aa(1:9, 1:9)
844  call calinverse(9, xj)
845  ! --- [Kaa]-1 * fa
846  icdisp(1:9) = matmul( xj(:, :), vect_a(1:3*ndof) )
847  ! --- fd - [Kda] * [Kaa]-1 * fa
848  vect(1:nn*ndof) = vect(1:nn*ndof) - matmul( stiff_da(1:nn*ndof, 1:3*ndof), icdisp(1:3*ndof) )
849 
850  end subroutine tload_c3d8ic
851 
852 
853 end module m_static_lib_3dic
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
elementinfo::getshapefunc
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coordinate system.
Definition: element.f90:647
elementinfo::getjacobian
subroutine getjacobian(fetype, nn, localcoord, elecoord, det, jacobian, inverse)
calculate Jacobian matrix, its determinant and inverse
Definition: element.f90:820
elementinfo::getquadpoint
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:489
elementinfo::getglobalderiv
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:741
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
hecmw
Definition: hecmw.f90:6
m_common_struct
This modules defines common structures for fem analysis.
Definition: m_common_struct.f90:6
m_common_struct::g_localcoordsys
type(tlocalcoordsys), dimension(:), pointer, save g_localcoordsys
Definition: m_common_struct.f90:19
m_common_struct::set_localcoordsys
subroutine set_localcoordsys(coords, coordsys, outsys, ierr)
setup of coordinate system
Definition: m_common_struct.f90:72
m_elastoplastic
This module provide functions for elastoplastic calculation.
Definition: Elastoplastic.f90:6
m_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.f90:15
m_fstr::ref_temp
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.f90:136
m_matmatrix
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
m_matmatrix::matlmatrix
subroutine matlmatrix(gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp, hdflag)
Calculate constituive matrix.
Definition: calMatMatrix.f90:30
m_static_lib_3d
This module provides common functions of Solid elements.
Definition: static_LIB_3d.f90:6
m_static_lib_3d::geomat_c3
subroutine geomat_c3(stress, mat)
Definition: static_LIB_3d.f90:19
m_static_lib_3d::cal_thermal_expansion_c3
subroutine cal_thermal_expansion_c3(tt0, ttc, material, coordsys, matlaniso, EPSTH)
Definition: static_LIB_3d.f90:519
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_3dic
Eight-node hexagonal element with incompatible mode.
Definition: static_LIB_3dIC.f90:9
m_static_lib_3dic::stf_c3d8ic
subroutine stf_c3d8ic(etype, nn, ecoord, gausses, stiff, cdsys_ID, coords, time, tincr, u, aux, temperature)
CALCULATION STIFF Matrix for C3D8IC ELEMENT.
Definition: static_LIB_3dIC.f90:24
m_static_lib_3dic::update_c3d8ic
subroutine update_c3d8ic(etype, nn, ecoord, u, du, ddu, cdsys_ID, coords, qf, gausses, iter, time, tincr, aux, ddaux, TT, T0, TN)
Update strain and stress inside element.
Definition: static_LIB_3dIC.f90:223
m_static_lib_3dic::tload_c3d8ic
subroutine tload_c3d8ic(etype, nn, xx, yy, zz, tt, t0, gausses, VECT, cdsys_ID, coords)
This subroutine calculates thermal loading.
Definition: static_LIB_3dIC.f90:681
m_utilities
This module provides aux functions.
Definition: utilities.f90:6
m_utilities::transformation
subroutine transformation(jacob, tm)
Definition: utilities.f90:376
m_utilities::calinverse
subroutine calinverse(NN, A)
calculate inverse of matrix a
Definition: utilities.f90:295
mmaterial
This module summarizes all information of material properties.
Definition: material.f90:6
mmaterial::d3
integer(kind=kint), parameter d3
Definition: material.f90:84
mmaterial::totallag
integer(kind=kint), parameter totallag
Definition: material.f90:14
mmaterial::mc_orthoexp
character(len=dict_key_length) mc_orthoexp
Definition: material.f90:144
mmaterial::infinitesimal
integer(kind=kint), parameter infinitesimal
Definition: material.f90:13
mmaterial::updatelag
integer(kind=kint), parameter updatelag
Definition: material.f90:15
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
mmechgauss::tgaussstatus
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13