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