FrontISTR  5.2.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 !-------------------------------------------------------------------------------
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), optional :: 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 imcompatible 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  if( present(temperature) ) then
102  call getshapefunc( etype, naturalcoord, spfunc )
103  temp = dot_product( temperature, spfunc )
104  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
105  else
106  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys )
107  end if
108 
109  if( flag == updatelag ) then
110  call geomat_c3( gausses(lx)%stress, mat )
111  d(:, :) = d(:, :)-mat
112  endif
113 
114  ! -- Derivative of shape function of imcompatible mode --
115  ! [ -2*a 0, 0 ]
116  ! [ 0, -2*b, 0 ]
117  ! [ 0, 0, -2*c ]
118  ! we don't call getShapeDeriv but use above value directly for computation efficiency
119  gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse(1, :)/det
120  gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse(2, :)/det
121  gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse(3, :)/det
122 
123  wg = getweight(etype, lx)*det
124  b(1:6, 1:(nn+3)*ndof)=0.0d0
125  do j = 1, nn+3
126  b(1, 3*j-2) = gderiv(j, 1)
127  b(2, 3*j-1) = gderiv(j, 2)
128  b(3, 3*j ) = gderiv(j, 3)
129  b(4, 3*j-2) = gderiv(j, 2)
130  b(4, 3*j-1) = gderiv(j, 1)
131  b(5, 3*j-1) = gderiv(j, 3)
132  b(5, 3*j ) = gderiv(j, 2)
133  b(6, 3*j-2) = gderiv(j, 3)
134  b(6, 3*j ) = gderiv(j, 1)
135  end do
136 
137  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
138  if( flag == totallag ) then
139  ! ---dudx(i,j) ==> gdispderiv(i,j)
140  gdispderiv(1:ndof,1:ndof) = matmul( unode(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) )
141  do j = 1, nn+3
142 
143  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
144  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
145  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
146  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
147  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
148  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
149  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
150  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
151  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
152  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
153  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
154  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
155  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
156  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
157  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
158  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
159  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
160  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
161 
162  end do
163  ! ---BL = BL0 + BL1
164  do j = 1, (nn+3)*ndof
165  b(:, j) = b(:, j)+b1(:, j)
166  end do
167 
168  end if
169 
170  db(1:6, 1:(nn+3)*ndof) = matmul( d, b(1:6, 1:(nn+3)*ndof) )
171  forall( i=1:(nn+3)*ndof, j=1:(nn+3)*ndof )
172  tmpstiff(i, j) = tmpstiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
173  end forall
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  forall( i=1:(nn+3)*ndof, j=1:(nn+3)*ndof )
203  tmpstiff(i, j) = tmpstiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
204  end forall
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), optional :: tt(nn)
249  real(kind=kreal), intent(in), optional :: t0(nn)
250  real(kind=kreal), intent(in), optional :: tn(nn)
251 
252  ! LCOAL VARIAVLES
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 imcompatible 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  if( present(tt) ) then
305  call getshapefunc( etype, naturalcoord, spfunc )
306  ttc = dot_product( tt, spfunc )
307  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, ttc )
308  else
309  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys )
310  end if
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 imcompatible 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  forall( i=1:3*ndof, j=1:nn*ndof )
375  stiff_ad(i, j) = stiff_ad(i, j)+dot_product( b(:, nn*ndof+i), db(:, j) )*wg
376  end forall
377  forall( i=1:3*ndof, j=1:3*ndof )
378  stiff_aa(i, j) = stiff_aa(i, j)+dot_product( b(:, nn*ndof+i), db(:, nn*ndof+j) )*wg
379  end forall
380 
381  ! calculate the stress matrix ( TOTAL LAGRANGE METHOD )
382  if( flag == totallag .OR. flag == updatelag ) then
383  stress(1:6) = gausses(lx)%stress
384  do j = 1, nn+3
385  bn(1, 3*j-2) = gderiv(j, 1)
386  bn(2, 3*j-1) = gderiv(j, 1)
387  bn(3, 3*j ) = gderiv(j, 1)
388  bn(4, 3*j-2) = gderiv(j, 2)
389  bn(5, 3*j-1) = gderiv(j, 2)
390  bn(6, 3*j ) = gderiv(j, 2)
391  bn(7, 3*j-2) = gderiv(j, 3)
392  bn(8, 3*j-1) = gderiv(j, 3)
393  bn(9, 3*j ) = gderiv(j, 3)
394  end do
395  smat(:, :) = 0.0d0
396  do j = 1, 3
397  smat(j , j ) = stress(1)
398  smat(j , j+3) = stress(4)
399  smat(j , j+6) = stress(6)
400  smat(j+3, j ) = stress(4)
401  smat(j+3, j+3) = stress(2)
402  smat(j+3, j+6) = stress(5)
403  smat(j+6, j ) = stress(6)
404  smat(j+6, j+3) = stress(5)
405  smat(j+6, j+6) = stress(3)
406  end do
407  sbn(1:9, 1:(nn+3)*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:(nn+3)*ndof) )
408  forall( i=1:3*ndof, j=1:nn*ndof )
409  stiff_ad(i, j) = stiff_ad(i, j)+dot_product( bn(:, nn*ndof+i), sbn(:, j) )*wg
410  end forall
411  forall( i=1:3*ndof, j=1:3*ndof )
412  stiff_aa(i, j) = stiff_aa(i, j)+dot_product( bn(:, nn*ndof+i), sbn(:, nn*ndof+j) )*wg
413  end forall
414  end if
415 
416  end do
417 
418  ! -----recover incompatible dof of nodal displacement
419  xj(1:3*ndof, 1:3*ndof)= stiff_aa(1:3*ndof, 1:3*ndof)
420  call calinverse(3*ndof, xj)
421  ! --- [Kad] * ddu
422  do i=1,nn
423  tmp_d((i-1)*ndof+1:i*ndof) = ddu(1:ndof, i)
424  enddo
425  tmpforce(1:3*ndof) = matmul( stiff_ad(1:3*ndof, 1:nn*ndof), tmp_d(1:nn*ndof) )
426  ! --- ddaux = -[Kaa]-1 * ([Kad] * ddu)
427  tmp_a(1:3*ndof) = -matmul( xj(1:3*ndof, 1:3*ndof), tmpforce(1:3*ndof) )
428  do i=1,3
429  ddaux(1:ndof, i) = tmp_a((i-1)*ndof+1:i*ndof)
430  enddo
431 
432 
433  !---------------------------------------------------------------------
434 
435  totaldisp(:,1:nn) = u(:,:)+ du(:,:)
436  totaldisp(:,nn+1:nn+3) = aux(:,:) + ddaux(:,:)
437 
438  !---------------------------------------------------------------------
439 
440  qf(:) = 0.0d0
441  qf_a(:) = 0.0d0
442 
443  stiff_ad(:, :) = 0.0d0
444  stiff_aa(:, :) = 0.0d0
445 
446  !---------------------------------------------------------------------
447 
448  if( flag == updatelag ) then
449  elem(:,:) = (0.5d0*du(:,:)+u(:,:) ) +ecoord(:,:)
450  elem1(:,:) = (du(:,:)+u(:,:) ) +ecoord(:,:)
451  ! elem = elem1
452  totaldisp(:,1:nn) = du(:,:)
453  if( iter == 1 ) aux(:,:) = 0.0d0 !--- is this correct???
454  totaldisp(:,nn+1:nn+3) = aux(:,:)
455  end if
456 
457  matlaniso = .false.
458  if( present(tt) .and. cdsys_id > 0 ) then
459  ina = tt(1)
460  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
461  if( .not. ierr ) matlaniso = .true.
462  end if
463 
464  ! --- Inverse of Jacobian at elemental center
465  naturalcoord(:) = 0.0d0
466  call getjacobian(etype, nn, naturalcoord, elem, det, jacobian, inverse)
467  inverse(:, :)= inverse(:, :)*det
468  if( flag == updatelag ) then
469  call getjacobian(etype, nn, naturalcoord, elem1, det, jacobian, inverse1)
470  inverse1(:, :)= inverse1(:, :)*det
471  call getjacobian(etype, nn, naturalcoord, ecoord, det, jacobian, inverse0)
472  inverse0(:, :)= inverse0(:, :)*det !for deformation gradient F
473  endif
474 
475  do lx = 1, numofquadpoints(etype)
476 
477  mtype = gausses(lx)%pMaterial%mtype
478 
479  call getquadpoint( etype, lx, naturalcoord(:) )
480  call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
481 
482  if( cdsys_id > 0 ) then
483  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:,:), serr )
484  if( serr == -1 ) stop "Fail to setup local coordinate"
485  if( serr == -2 ) then
486  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
487  end if
488  end if
489 
490  ! ========================================================
491  ! UPDATE STRAIN and STRESS
492  ! ========================================================
493 
494  ! Thermal Strain
495  epsth = 0.0d0
496  if( present(tt) .AND. present(t0) ) then
497  call getshapefunc(etype, naturalcoord, spfunc)
498  ttc = dot_product(tt, spfunc)
499  tt0 = dot_product(t0, spfunc)
500  ttn = dot_product(tn, spfunc)
501  call cal_thermal_expansion_c3( tt0, ttc, gausses(lx)%pMaterial, coordsys, matlaniso, epsth )
502  end if
503 
504  ! -- Derivative of shape function of imcompatible mode --
505  ! [ -2*a 0, 0 ]
506  ! [ 0, -2*b, 0 ]
507  ! [ 0, 0, -2*c ]
508  ! we don't call getShapeDeriv but use above value directly for computation efficiency
509  gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse(1, :)/det
510  gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse(2, :)/det
511  gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse(3, :)/det
512 
513  ! Small strain
514  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) )
515  dstrain(1) = gdispderiv(1, 1)
516  dstrain(2) = gdispderiv(2, 2)
517  dstrain(3) = gdispderiv(3, 3)
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(:) ! allright?
522 
523  f(1:3,1:3) = 0.d0; f(1,1)=1.d0; f(2,2)=1.d0; f(3,3)=1.d0; !deformation gradient
524  if( flag == infinitesimal ) then
525  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
526  f(1:3,1:3) = f(1:3,1:3) + gdispderiv(1:3,1:3)
527 
528  else if( flag == totallag ) then
529  ! Green-Lagrange strain
530  dstrain(1) = dstrain(1)+0.5d0*dot_product( gdispderiv(:, 1), gdispderiv(:, 1) )
531  dstrain(2) = dstrain(2)+0.5d0*dot_product( gdispderiv(:, 2), gdispderiv(:, 2) )
532  dstrain(3) = dstrain(3)+0.5d0*dot_product( gdispderiv(:, 3), gdispderiv(:, 3) )
533  dstrain(4) = dstrain(4)+( gdispderiv(1, 1)*gdispderiv(1, 2) &
534  +gdispderiv(2, 1)*gdispderiv(2, 2)+gdispderiv(3, 1)*gdispderiv(3, 2) )
535  dstrain(5) = dstrain(5)+( gdispderiv(1, 2)*gdispderiv(1, 3) &
536  +gdispderiv(2, 2)*gdispderiv(2, 3)+gdispderiv(3, 2)*gdispderiv(3, 3) )
537  dstrain(6) = dstrain(6)+( gdispderiv(1, 1)*gdispderiv(1, 3) &
538  +gdispderiv(2, 1)*gdispderiv(2, 3)+gdispderiv(3, 1)*gdispderiv(3, 3) )
539 
540  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
541 
542  else if( flag == updatelag ) then
543  ! CALL GEOMAT_C3( gausses(LX)%stress, mat )
544  ! D(:, :) = D(:, :)+mat(:, :)
545  rot = 0.0d0
546  rot(1, 2)= 0.5d0*(gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
547  rot(2, 3)= 0.5d0*(gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
548  rot(1, 3)= 0.5d0*(gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
549 
550  gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+ dstrain(1:6)+epsth(:)
551  call getglobalderiv(etype, nn, naturalcoord, ecoord, det0, gderiv0)
552  gderiv0(nn+1, :) = -2.0d0*naturalcoord(1)*inverse0(1, :)/det0
553  gderiv0(nn+2, :) = -2.0d0*naturalcoord(2)*inverse0(2, :)/det0
554  gderiv0(nn+3, :) = -2.0d0*naturalcoord(3)*inverse0(3, :)/det0
555  f(1:3,1:3) = f(1:3,1:3) + matmul( u(1:ndof, 1:nn+3)+du(1:ndof, 1:nn+3), gderiv0(1:nn+3, 1:ndof) )
556 
557  end if
558 
559  ! Update stress
560  if( present(tt) .AND. present(t0) ) then
561  call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn )
562  else
563  call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr )
564  end if
565 
566  ! ========================================================
567  ! calculate the internal force ( equivalent nodal force )
568  ! ========================================================
569  ! Small strain
570  b(1:6, 1:(nn+3)*ndof) = 0.0d0
571  do j=1,nn+3
572  b(1,3*j-2) = gderiv(j, 1)
573  b(2,3*j-1) = gderiv(j, 2)
574  b(3,3*j ) = gderiv(j, 3)
575  b(4,3*j-2) = gderiv(j, 2)
576  b(4,3*j-1) = gderiv(j, 1)
577  b(5,3*j-1) = gderiv(j, 3)
578  b(5,3*j ) = gderiv(j, 2)
579  b(6,3*j-2) = gderiv(j, 3)
580  b(6,3*j ) = gderiv(j, 1)
581  end do
582 
583  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
584  if( flag == infinitesimal ) then
585 
586  else if( flag == totallag ) then
587 
588  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn+3), gderiv(1:nn+3, 1:ndof) )
589  b1(1:6, 1:(nn+3)*ndof)=0.0d0
590  do j = 1,nn+3
591  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
592  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
593  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
594  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
595  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
596  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
597  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
598  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
599  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
600  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
601  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
602  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
603  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
604  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
605  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
606  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
607  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
608  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
609  end do
610  ! BL = BL0 + BL1
611  do j=1,(nn+3)*ndof
612  b(:,j) = b(:,j)+b1(:,j)
613  end do
614 
615  else if( flag == updatelag ) then
616 
617  call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv(1:nn,1:3))
618 
619  ! -- Derivative of shape function of imcompatible mode --
620  ! [ -2*a 0, 0 ]
621  ! [ 0, -2*b, 0 ]
622  ! [ 0, 0, -2*c ]
623  ! we don't call getShapeDeriv but use above value directly for computation efficiency
624  gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse1(1, :)/det
625  gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse1(2, :)/det
626  gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse1(3, :)/det
627 
628  b(1:6, 1:(nn+3)*ndof) = 0.0d0
629  do j = 1, nn+3
630  b(1, 3*j-2) = gderiv(j, 1)
631  b(2, 3*j-1) = gderiv(j, 2)
632  b(3, 3*j ) = gderiv(j, 3)
633  b(4, 3*j-2) = gderiv(j, 2)
634  b(4, 3*j-1) = gderiv(j, 1)
635  b(5, 3*j-1) = gderiv(j, 3)
636  b(5, 3*j ) = gderiv(j, 2)
637  b(6, 3*j-2) = gderiv(j, 3)
638  b(6, 3*j ) = gderiv(j, 1)
639  end do
640 
641  end if
642 
643  wg=getweight( etype, lx )*det
644 
645  db(1:6, 1:(nn+3)*ndof) = matmul( d, b(1:6, 1:(nn+3)*ndof) )
646  forall( i=1:3*ndof, j=1:nn*ndof )
647  stiff_ad(i, j) = stiff_ad(i, j)+dot_product( b(:, nn*ndof+i), db(:, j) )*wg
648  end forall
649  forall( i=1:3*ndof, j=1:3*ndof )
650  stiff_aa(i, j) = stiff_aa(i, j)+dot_product( b(:, nn*ndof+i), db(:, nn*ndof+j) )*wg
651  end forall
652 
653  ! calculate the Internal Force
654  qf(1:nn*ndof) &
655  = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
656  qf_a(1:3*ndof) &
657  = qf_a(1:3*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,nn*ndof+1:(nn+3)*ndof) )*wg
658  end do
659 
660  ! condence ( qf - [Kda] * [Kaa]-1 * qf_a )
661  xj(1:9, 1:9)= stiff_aa(1:3*ndof, 1:3*ndof)
662  call calinverse(9, xj)
663  tmpdisp(:) = matmul( xj(:,:), qf_a(:) )
664  do i=1,nn*ndof
665  qf(i) = qf(i) - dot_product( stiff_ad(:,i), tmpdisp(:) )
666  enddo
667  end subroutine update_c3d8ic
668 
669 
671  !----------------------------------------------------------------------*
672  subroutine tload_c3d8ic &
673  (etype, nn, xx, yy, zz, tt, t0, &
674  gausses, vect, cdsys_id, coords)
675  !----------------------------------------------------------------------*
676 
677  use m_fstr
678  use mmechgauss
679  use m_matmatrix
680  use m_common_struct
681 
682  !---------------------------------------------------------------------
683 
684  integer(kind=kint), parameter :: ndof=3
685  integer(kind=kint), intent(in) :: etype
686  integer(kind=kint), intent(in) :: nn
687  real(kind=kreal), intent(in) :: xx(nn), yy(nn), zz(nn)
688  real(kind=kreal), intent(in) :: tt(nn), t0(nn)
689  type(tgaussstatus), intent(inout) :: gausses(:)
690  real(kind=kreal), intent(out) :: vect(nn*ndof)
691  integer(kind=kint), intent(in) :: cdsys_id
692  real(kind=kreal), intent(inout) :: coords(3, 3)
693 
694  !---------------------------------------------------------------------
695 
696  real(kind=kreal) :: alp, alp0
697  real(kind=kreal) :: d(6, 6), b(6, ndof*(nn+3)), db_a(6, ndof*3)
698  real(kind=kreal) :: det, wg, ecoord(3, nn)
699  integer(kind=kint) :: i, j, ic, serr
700  real(kind=kreal) :: epsth(6), sgm(6), h(nn), alpo(3), alpo0(3), coordsys(3, 3), xj(9, 9)
701  real(kind=kreal) :: naturalcoord(3), gderiv(nn+3, 3)
702  real(kind=kreal) :: jacobian(3, 3),inverse(3, 3)
703  real(kind=kreal) :: stiff_da(nn*3, 3*3), stiff_aa(3*3, 3*3)
704  real(kind=kreal) :: vect_a(3*3)
705  real(kind=kreal) :: icdisp(9)
706  real(kind=kreal) :: tempc, temp0, thermal_eps, tm(6, 6), outa(1), ina(1)
707  logical :: ierr, matlaniso
708 
709  !---------------------------------------------------------------------
710 
711  matlaniso = .false.
712  if( cdsys_id > 0 ) then ! cannot define aniso exapansion when no local coord defined
713  ina = tt(1)
714  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
715  if( .not. ierr ) matlaniso = .true.
716  end if
717 
718  ecoord(1, :) = xx(:)
719  ecoord(2, :) = yy(:)
720  ecoord(3, :) = zz(:)
721  ! ---- Calculate enhanced displacement at first
722  ! --- Inverse of Jacobian at elemental center
723  naturalcoord(:) = 0.0d0
724  call getjacobian(etype, nn, naturalcoord, ecoord, det, jacobian, inverse)
725  inverse(:, :) = inverse(:, :)*det
726  ! ---- We now calculate stiff matrix include imcompatible mode
727  ! [ Kdd Kda ]
728  ! [ Kad Kaa ]
729  stiff_da(:, :) = 0.0d0
730  stiff_aa(:, :) = 0.0d0
731  b(1:6, 1:(nn+3)*ndof) = 0.0d0
732  vect(1:nn*ndof) = 0.0d0
733  vect_a(1:3*ndof) = 0.0d0
734 
735  do ic = 1, numofquadpoints(etype)
736 
737  call getquadpoint(etype, ic, naturalcoord)
738  call getglobalderiv( etype, nn, naturalcoord, ecoord, det, gderiv(1:nn, 1:3) )
739 
740  if( cdsys_id > 0 ) then
741  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
742  if( serr == -1 ) stop "Fail to setup local coordinate"
743  if( serr == -2 ) then
744  write(*,*) "WARNING! Cannot setup local coordinate, it is modified automatically"
745  end if
746  end if
747 
748  call getshapefunc( etype, naturalcoord, h(1:nn) )
749  tempc = dot_product( h(1:nn), tt(1:nn) )
750  temp0 = dot_product( h(1:nn), t0(1:nn) )
751  call matlmatrix( gausses(ic), d3, d, 1.d0, 1.0d0, coordsys, tempc )
752 
753  ! -- Derivative of shape function of imcompatible mode --
754  ! [ -2*a 0, 0 ]
755  ! [ 0, -2*b, 0 ]
756  ! [ 0, 0, -2*c ]
757  ! we don't call getShapeDeriv but use above value directly for computation efficiency
758  gderiv(nn+1, :) = -2.0d0*naturalcoord(1)*inverse(1, :)/det
759  gderiv(nn+2, :) = -2.0d0*naturalcoord(2)*inverse(2, :)/det
760  gderiv(nn+3, :) = -2.0d0*naturalcoord(3)*inverse(3, :)/det
761 
762  wg = getweight(etype, ic)*det
763  do j = 1, nn+3
764  b(1, 3*j-2) = gderiv(j, 1)
765  b(2, 3*j-1) = gderiv(j, 2)
766  b(3, 3*j ) = gderiv(j, 3)
767  b(4, 3*j-2) = gderiv(j, 2)
768  b(4, 3*j-1) = gderiv(j, 1)
769  b(5, 3*j-1) = gderiv(j, 3)
770  b(5, 3*j ) = gderiv(j, 2)
771  b(6, 3*j-2) = gderiv(j, 3)
772  b(6, 3*j ) = gderiv(j, 1)
773  end do
774 
775  db_a(1:6, 1:3*ndof) = matmul( d, b(1:6, nn*ndof+1:(nn+3)*ndof) )
776  forall( i=1:nn*ndof, j=1:3*ndof )
777  stiff_da(i, j) = stiff_da(i, j) + dot_product( b(1:6, i), db_a(1:6, j) ) * wg
778  end forall
779  forall( i=1:3*ndof, j=1:3*ndof )
780  stiff_aa(i, j) = stiff_aa(i, j) + dot_product( b(1:6, nn*ndof+i), db_a(1:6, j) ) * wg
781  end forall
782 
783  ina(1) = tempc
784  if( matlaniso ) then
785  call fetch_tabledata( mc_orthoexp, gausses(ic)%pMaterial%dict, alpo(:), ierr, ina )
786  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
787  else
788  call fetch_tabledata( mc_themoexp, gausses(ic)%pMaterial%dict, outa(:), ierr, ina )
789  if( ierr ) outa(1) = gausses(ic)%pMaterial%variables(m_exapnsion)
790  alp = outa(1)
791  end if
792  ina(1) = temp0
793  if( matlaniso ) then
794  call fetch_tabledata( mc_orthoexp, gausses(ic)%pMaterial%dict, alpo0(:), ierr, ina )
795  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
796  else
797  call fetch_tabledata( mc_themoexp, gausses(ic)%pMaterial%dict, outa(:), ierr, ina )
798  if( ierr ) outa(1) = gausses(ic)%pMaterial%variables(m_exapnsion)
799  alp0 = outa(1)
800  end if
801 
802  !**
803  !** THERMAL strain
804  !**
805  if( matlaniso ) then
806  do j = 1, 3
807  epsth(j) = alpo(j)*(tempc-ref_temp)-alpo0(j)*(temp0-ref_temp)
808  end do
809  epsth(4:6) = 0.0d0
810  call transformation(coordsys, tm)
811  epsth(:) = matmul( epsth(:), tm ) ! to global coord
812  epsth(4:6) = epsth(4:6)*2.0d0
813  else
814  thermal_eps = alp*(tempc-ref_temp)-alp0*(temp0-ref_temp)
815  epsth(1:3) = thermal_eps
816  epsth(4:6) = 0.0d0
817  end if
818 
819  !**
820  !** SET SGM {s}=[D]{e}
821  !**
822  sgm(:) = matmul( d(:, :), epsth(:) )
823 
824  !**
825  !** CALCULATE LOAD {F}=[B]T{e}
826  !**
827  vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(1:6), b(1:6, 1:nn*ndof) )*wg
828  vect_a(1:3*ndof) = vect_a(1:3*ndof)+matmul( sgm(1:6), b(1:6, nn*ndof+1:(nn+3)*ndof) )*wg
829 
830  end do
831 
832  ! --- condense vect
833  xj(1:9,1:9)= stiff_aa(1:9, 1:9)
834  call calinverse(9, xj)
835  ! --- [Kaa]-1 * fa
836  icdisp(1:9) = matmul( xj(:, :), vect_a(1:3*ndof) )
837  ! --- fd - [Kda] * [Kaa]-1 * fa
838  vect(1:nn*ndof) = vect(1:nn*ndof) - matmul( stiff_da(1:nn*ndof, 1:3*ndof), icdisp(1:3*ndof) )
839 
840  end subroutine tload_c3d8ic
841 
842 
843 end module m_static_lib_3dic
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
Definition: element.f90:640
subroutine getjacobian(fetype, nn, localcoord, elecoord, det, jacobian, inverse)
calculate Jacobian matrix, its determinant and inverse
Definition: element.f90:813
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:482
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:734
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:528
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:445
Definition: hecmw.f90:6
This modules defines common structures for fem analysis.
type(tlocalcoordsys), dimension(:), pointer, save g_localcoordsys
subroutine set_localcoordsys(coords, coordsys, outsys, ierr)
setup of coordinate system
This module provide functions for elastoplastic calculation.
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.f90:120
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
subroutine matlmatrix(gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp)
Calculate constituive matrix.
This module provide common functions of Solid elements.
subroutine geomat_c3(stress, mat)
subroutine cal_thermal_expansion_c3(tt0, ttc, material, coordsys, matlaniso, EPSTH)
subroutine update_stress3d(flag, gauss, rot, dstrain, F, coordsys, time, tincr, ttc, tt0, ttn)
Eight-node hexagonal element with imcompatible mode.
subroutine stf_c3d8ic(etype, nn, ecoord, gausses, stiff, cdsys_ID, coords, time, tincr, u, aux, temperature)
CALCULATION STIFF Matrix for C3D8IC ELEMENT.
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.
subroutine tload_c3d8ic(etype, nn, xx, yy, zz, tt, t0, gausses, VECT, cdsys_ID, coords)
This subroutine calculatess thermal loading.
This module provides aux functions.
Definition: utilities.f90:6
subroutine transformation(jacob, tm)
Definition: utilities.f90:339
subroutine calinverse(NN, A)
calculate inverse of matrix a
Definition: utilities.f90:258
This module summarizes all infomation of material properties.
Definition: material.f90:6
integer(kind=kint), parameter d3
Definition: material.f90:76
integer(kind=kint), parameter totallag
Definition: material.f90:14
character(len=dict_key_length) mc_orthoexp
Definition: material.f90:123
integer(kind=kint), parameter infinitesimal
Definition: material.f90:13
integer(kind=kint), parameter updatelag
Definition: material.f90:15
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13