FrontISTR  5.9.0
Large-scale structural analysis program with finit element method
static_LIB_3d_up.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 !-------------------------------------------------------------------------------
6 module m_static_lib_3d_up
7 
8  use hecmw, only : kint, kreal
9  use elementinfo
10 
11  implicit none
12 
13 contains
14 
16  subroutine stf_c3_up &
17  (etype, nn, ecoord, gausses, stiff, cdsys_id, coords, &
18  time, tincr, nn_p, lambda, u, temperature)
19 
20  use mmechgauss
21  use m_matmatrix
22  use m_common_struct
23  use m_static_lib_3d, only: geomat_c3
24 
25  integer(kind=kint), intent(in) :: etype
26  integer(kind=kint), intent(in) :: nn
27  real(kind=kreal), intent(in) :: ecoord(3, nn)
28  type(tgaussstatus), intent(in) :: gausses(:)
29  real(kind=kreal), intent(out) :: stiff(:, :)
30  integer(kind=kint), intent(in) :: cdsys_id
31  real(kind=kreal), intent(inout) :: coords(3, 3)
32  real(kind=kreal), intent(in) :: time
33  real(kind=kreal), intent(in) :: tincr
34  integer(kind=kint), intent(in) :: nn_p
35  real(kind=kreal), intent(in) :: lambda(nn_p)
36  real(kind=kreal), intent(in), optional :: u(:, :)
37  real(kind=kreal), intent(in), optional :: temperature(nn)
38 
39  integer(kind=kint) :: flag
40  integer(kind=kint), parameter :: ndof = 3
41  real(kind=kreal) :: d(6, 6), b(6, ndof*nn), db(6, ndof*nn)
42  real(kind=kreal) :: gderiv(nn, 3), stress(6), mat(6, 6)
43  real(kind=kreal) :: det, wg, temp, spfunc(nn)
44  integer(kind=kint) :: i, j, lx, serr
45  real(kind=kreal) :: naturalcoord(3), coordsys(3, 3)
46  real(kind=kreal) :: gdispderiv(3, 3)
47  real(kind=kreal) :: b1(6, ndof*nn)
48  real(kind=kreal) :: smat(9, 9), elem(3, nn)
49  real(kind=kreal) :: bn(9, ndof*nn), sbn(9, ndof*nn)
50  integer(kind=kint) :: na, nb
51  integer(kind=kint) :: na_p, nb_p
52  integer(kind=kint) :: isize, jsize
53  integer(kind=kint) :: jsize1, jsize2, jsize3
54  real(kind=kreal) :: stiff_up(3*nn, nn_p)
55  real(kind=kreal) :: stiff_pp(nn_p, nn_p), stiff_pp_inv(nn_p, nn_p)
56  real(kind=kreal) :: stiff_up_stiff_pp_inv(3*nn, nn_p)
57  real(kind=kreal) :: stiff_up_stiff_pp_inv_stiff_up(3*nn, 3*nn)
58  real(kind=kreal) :: alpha_inv
59  real(kind=kreal) :: g
60  real(kind=kreal) :: d2(6)
61  real(kind=kreal) :: bd2(3*nn)
62  real(kind=kreal) :: lambda_lx
63  real(kind=kreal) :: spfunc_p(nn_p)
64 
65  stiff(:, :) = 0.0d0
66  stiff_up(:, :) = 0.0d0
67  stiff_pp(:, :) = 0.0d0
68 
69  ! we suppose the same material type in the element
70  flag = gausses(1)%pMaterial%nlgeom_flag
71  if( .not. present(u) ) flag = infinitesimal ! enforce to infinitesimal deformation analysis
72  elem(:, :) = ecoord(:, :)
73  ! Updated Lagrangian : evaluate on the current (deformed) configuration
74  if( flag == updatelag ) elem(:, :) = ecoord(:, :)+u(:, :)
75  if( flag == infinitesimal ) gdispderiv(:, :) = 0.0d0
76 
77  coordsys = 0.0d0
78  coordsys(1, 1) = 1.0d0; coordsys(2, 2) = 1.0d0; coordsys(3, 3) = 1.0d0
79 
80  do lx = 1, numofquadpoints(etype)
81 
82  call getquadpoint( etype, lx, naturalcoord(:) )
83  call getglobalderiv( etype, nn, naturalcoord, elem, det, gderiv )
84 
85  if( cdsys_id > 0 ) then
86  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
87  if( serr == -1 ) stop "Fail to setup local coordinate"
88  if( serr == -2 ) write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
89  end if
90 
91  ! Displacement gradient (B1 nonlinear term : Total Lagrangian only)
92  if( flag == totallag ) then
93  gdispderiv(1:3, 1:3) = matmul( u(1:3, 1:nn), gderiv(1:nn, 1:3) )
94  else
95  gdispderiv(1:3, 1:3) = 0.0d0
96  end if
97 
98  if( nn_p == 1 ) spfunc_p(1) = 1.0d0
99 
100  ! Lagrange multiplier
101  lambda_lx = 0.0d0
102  do na_p = 1, nn_p
103  lambda_lx = lambda_lx+spfunc_p(na_p)*lambda(na_p)
104  end do
105 
106  if( present( temperature ) ) then
107  call getshapefunc( etype, naturalcoord, spfunc )
108  temp = dot_product( temperature, spfunc )
109  call matlmatrix_up( gausses(lx), d3, d, lambda_lx, d2, alpha_inv, g, time, tincr, coordsys, temp )
110  else
111  call matlmatrix_up( gausses(lx), d3, d, lambda_lx, d2, alpha_inv, g, time, tincr, coordsys )
112  end if
113 
114  ! Updated Lagrangian : material tangent correction with the Cauchy stress
115  if( flag == updatelag ) then
116  call geomat_c3( gausses(lx)%stress, mat )
117  d(:, :) = d(:, :)-mat
118  end if
119 
120  wg = getweight( etype, lx )*det
121 
122  do nb = 1, nn
123  jsize1 = 3*(nb-1)+1
124  jsize2 = 3*(nb-1)+2
125  jsize3 = 3*(nb-1)+3
126  b(1, jsize1) = gderiv(nb, 1)
127  b(2, jsize1) = 0.0d0
128  b(3, jsize1) = 0.0d0
129  b(4, jsize1) = gderiv(nb, 2)
130  b(5, jsize1) = 0.0d0
131  b(6, jsize1) = gderiv(nb, 3)
132  b(1, jsize2) = 0.0d0
133  b(2, jsize2) = gderiv(nb, 2)
134  b(3, jsize2) = 0.0d0
135  b(4, jsize2) = gderiv(nb, 1)
136  b(5, jsize2) = gderiv(nb, 3)
137  b(6, jsize2) = 0.0d0
138  b(1, jsize3) = 0.0d0
139  b(2, jsize3) = 0.0d0
140  b(3, jsize3) = gderiv(nb, 3)
141  b(4, jsize3) = 0.0d0
142  b(5, jsize3) = gderiv(nb, 2)
143  b(6, jsize3) = gderiv(nb, 1)
144  end do
145 
146  do nb = 1, nn
147  jsize1 = 3*(nb-1)+1
148  jsize2 = 3*(nb-1)+2
149  jsize3 = 3*(nb-1)+3
150  b1(1, jsize1) = gdispderiv(1, 1)*gderiv(nb, 1)
151  b1(2, jsize1) = gdispderiv(1, 2)*gderiv(nb, 2)
152  b1(3, jsize1) = gdispderiv(1, 3)*gderiv(nb, 3)
153  b1(4, jsize1) = gdispderiv(1, 2)*gderiv(nb, 1)+gdispderiv(1, 1)*gderiv(nb, 2)
154  b1(5, jsize1) = gdispderiv(1, 2)*gderiv(nb, 3)+gdispderiv(1, 3)*gderiv(nb, 2)
155  b1(6, jsize1) = gdispderiv(1, 3)*gderiv(nb, 1)+gdispderiv(1, 1)*gderiv(nb, 3)
156  b1(1, jsize2) = gdispderiv(2, 1)*gderiv(nb, 1)
157  b1(2, jsize2) = gdispderiv(2, 2)*gderiv(nb, 2)
158  b1(3, jsize2) = gdispderiv(2, 3)*gderiv(nb, 3)
159  b1(4, jsize2) = gdispderiv(2, 2)*gderiv(nb, 1)+gdispderiv(2, 1)*gderiv(nb, 2)
160  b1(5, jsize2) = gdispderiv(2, 2)*gderiv(nb, 3)+gdispderiv(2, 3)*gderiv(nb, 2)
161  b1(6, jsize2) = gdispderiv(2, 3)*gderiv(nb, 1)+gdispderiv(2, 1)*gderiv(nb, 3)
162  b1(1, jsize3) = gdispderiv(3, 1)*gderiv(nb, 1)
163  b1(2, jsize3) = gdispderiv(3, 2)*gderiv(nb, 2)
164  b1(3, jsize3) = gdispderiv(3, 3)*gderiv(nb, 3)
165  b1(4, jsize3) = gdispderiv(3, 2)*gderiv(nb, 1)+gdispderiv(3, 1)*gderiv(nb, 2)
166  b1(5, jsize3) = gdispderiv(3, 2)*gderiv(nb, 3)+gdispderiv(3, 3)*gderiv(nb, 2)
167  b1(6, jsize3) = gdispderiv(3, 3)*gderiv(nb, 1)+gdispderiv(3, 1)*gderiv(nb, 3)
168  end do
169 
170  do jsize = 1, 3*nn
171  b(:, jsize) = b(:, jsize)+b1(:, jsize)
172  end do
173 
174  db(1:6, 1:3*nn) = matmul( d(1:6, 1:6), b(1:6, 1:3*nn) )
175 
176  ! [ Kuu ] matrix
177  forall( isize = 1:3*nn, jsize = 1:3*nn )
178  stiff(isize, jsize) = stiff(isize, jsize)+wg*dot_product( b(:, isize), db(:, jsize) )
179  end forall
180 
181  ! calculate the initial stress matrix (Total/Updated Lagrangian)
182  if( flag == totallag .or. flag == updatelag ) then
183  stress(1:6) = gausses(lx)%stress
184 
185  do nb = 1, nn
186  jsize1 = 3*(nb-1)+1
187  jsize2 = 3*(nb-1)+2
188  jsize3 = 3*(nb-1)+3
189  bn(1, jsize1) = gderiv(nb, 1)
190  bn(2, jsize1) = 0.0d0
191  bn(3, jsize1) = 0.0d0
192  bn(4, jsize1) = gderiv(nb, 2)
193  bn(5, jsize1) = 0.0d0
194  bn(6, jsize1) = 0.0d0
195  bn(7, jsize1) = gderiv(nb, 3)
196  bn(8, jsize1) = 0.0d0
197  bn(9, jsize1) = 0.0d0
198  bn(1, jsize2) = 0.0d0
199  bn(2, jsize2) = gderiv(nb, 1)
200  bn(3, jsize2) = 0.0d0
201  bn(4, jsize2) = 0.0d0
202  bn(5, jsize2) = gderiv(nb, 2)
203  bn(6, jsize2) = 0.0d0
204  bn(7, jsize2) = 0.0d0
205  bn(8, jsize2) = gderiv(nb, 3)
206  bn(9, jsize2) = 0.0d0
207  bn(1, jsize3) = 0.0d0
208  bn(2, jsize3) = 0.0d0
209  bn(3, jsize3) = gderiv(nb, 1)
210  bn(4, jsize3) = 0.0d0
211  bn(5, jsize3) = 0.0d0
212  bn(6, jsize3) = gderiv(nb, 2)
213  bn(7, jsize3) = 0.0d0
214  bn(8, jsize3) = 0.0d0
215  bn(9, jsize3) = gderiv(nb, 3)
216  end do
217 
218  smat(:, :) = 0.0d0
219  do j = 1, 3
220  smat(j , j ) = stress(1)
221  smat(j , j+3) = stress(4)
222  smat(j , j+6) = stress(6)
223  smat(j+3, j ) = stress(4)
224  smat(j+3, j+3) = stress(2)
225  smat(j+3, j+6) = stress(5)
226  smat(j+6, j ) = stress(6)
227  smat(j+6, j+3) = stress(5)
228  smat(j+6, j+6) = stress(3)
229  end do
230 
231  sbn(1:9, 1:3*nn) = matmul( smat(1:9, 1:9), bn(1:9, 1:3*nn) )
232 
233  forall( isize = 1:3*nn, jsize = 1:3*nn )
234  stiff(isize, jsize) = stiff(isize, jsize)+wg*dot_product( bn(:, isize), sbn(:, jsize) )
235  end forall
236  end if
237 
238  do isize = 1, 3*nn
239  bd2(isize) = dot_product( b(:, isize), d2(:) )
240  end do
241 
242  ! [ Kup ] matrix
243  forall( isize = 1:3*nn, nb_p = 1:nn_p )
244  stiff_up(isize, nb_p) = stiff_up(isize, nb_p)+wg*bd2(isize)*spfunc_p(nb_p)
245  end forall
246 
247  ! [ Kpp ] matrix
248  forall( na_p = 1:nn_p, nb_p = 1:nn_p )
249  stiff_pp(na_p, nb_p) = stiff_pp(na_p, nb_p)-wg*alpha_inv*spfunc_p(na_p)*spfunc_p(nb_p)
250  end forall
251 
252  end do ! gauss loop
253 
254  ! [ Kpp ]^-1 matrix (constant pressure: nn_p = 1)
255  if( nn_p == 1 ) then
256  stiff_pp_inv(1, 1) = 1.0d0/stiff_pp(1, 1)
257  else
258  write(6, *) 'Error: nn_p should be equal to 1.'
259  return
260  end if
261 
262  stiff_up_stiff_pp_inv(1:3*nn, 1:nn_p) = matmul( stiff_up(1:3*nn, 1:nn_p), stiff_pp_inv(1:nn_p, 1:nn_p) )
263 
264  ! [ Kup ] [ Kpp ]^-1 [ Kup ]^T matrix
265  forall( isize = 1:3*nn, jsize = 1:3*nn )
266  stiff_up_stiff_pp_inv_stiff_up(isize, jsize) = dot_product( stiff_up_stiff_pp_inv(isize, :), stiff_up(jsize, :) )
267  end forall
268 
269  ! static condensation : [ Kuu ] - [ Kup ] [ Kpp ]^-1 [ Kup ]^T
270  forall( isize = 1:3*nn, jsize = 1:3*nn )
271  stiff(isize, jsize) = stiff(isize, jsize)-stiff_up_stiff_pp_inv_stiff_up(isize, jsize)
272  end forall
273 
274  end subroutine stf_c3_up
275 
276 
278  subroutine update_c3_up &
279  (etype, nn, ecoord, u, du, ddu, cdsys_id, coords, qf, &
280  gausses, iter, time, tincr, &
281  nn_p, lambda, ddlambda, tt, t0)
282 
283  use mmaterial
284  use mmechgauss
285  use m_matmatrix
286  use m_elastoplastic
287  use mhyperelastic
288  use m_common_struct
289  use m_utilities
290 
291  integer(kind=kint), intent(in) :: etype
292  integer(kind=kint), intent(in) :: nn
293  real(kind=kreal), intent(in) :: ecoord(3, nn)
294  real(kind=kreal), intent(in) :: u(3, nn)
295  real(kind=kreal), intent(in) :: du(3, nn)
296  real(kind=kreal), intent(in) :: ddu(3, nn)
297  integer(kind=kint), intent(in) :: cdsys_id
298  real(kind=kreal), intent(inout) :: coords(3, 3)
299  real(kind=kreal), intent(out) :: qf(nn*3)
300  type(tgaussstatus), intent(inout) :: gausses(:)
301  integer, intent(in) :: iter
302  real(kind=kreal), intent(in) :: time
303  real(kind=kreal), intent(in) :: tincr
304  integer(kind=kint), intent(in) :: nn_p
305  real(kind=kreal), intent(in) :: lambda(nn_p)
306  real(kind=kreal), intent(inout) :: ddlambda(nn_p)
307  real(kind=kreal), intent(in), optional :: tt(nn)
308  real(kind=kreal), intent(in), optional :: t0(nn)
309 
310  integer(kind=kint) :: flag
311  integer(kind=kint), parameter :: ndof = 3
312  real(kind=kreal) :: d(6, 6), b(6, ndof*nn), b1(6, ndof*nn)
313  real(kind=kreal) :: gderiv(nn, 3), gdispderiv(3, 3), det, wg
314  integer(kind=kint) :: i, j, lx, mtype, serr
315  real(kind=kreal) :: naturalcoord(3), rot(3, 3), spfunc(nn), coordsys(3, 3)
316  real(kind=kreal) :: totaldisp(3, nn), elem(3, nn), elem1(3, nn)
317  real(kind=kreal) :: dstrain(6), dstress(6), dumstress(3, 3), dum(3, 3)
318  real(kind=kreal) :: trd, p_bak
319  real(kind=kreal) :: ttc, tt0, outa(1), ina(1), epsth(6)
320  logical :: ierr
321  integer(kind=kint) :: na, nb
322  integer(kind=kint) :: na_p, nb_p
323  integer(kind=kint) :: isize, jsize
324  integer(kind=kint) :: jsize1, jsize2, jsize3
325  real(kind=kreal) :: totallambda(nn_p)
326  real(kind=kreal) :: stiff_up(3*nn, nn_p)
327  real(kind=kreal) :: stiff_pp(nn_p, nn_p), stiff_pp_inv(nn_p, nn_p)
328  real(kind=kreal) :: stiff_up_stiff_pp_inv(3*nn, nn_p)
329  real(kind=kreal) :: stiff_up_stiff_pp_inv_stiff_up(3*nn, 3*nn)
330  real(kind=kreal) :: stiff_up_stiff_pp_inv_qf_p(3*nn)
331  real(kind=kreal) :: alpha_inv
332  real(kind=kreal) :: g
333  real(kind=kreal) :: qf_p(nn_p)
334  real(kind=kreal) :: stiff_up_ddu(nn_p)
335  real(kind=kreal) :: d2(6)
336  real(kind=kreal) :: bd2(3*nn)
337  real(kind=kreal) :: lambda_lx
338  real(kind=kreal) :: spfunc_p(nn_p)
339 
340  stiff_up(:, :) = 0.0d0
341  stiff_pp(:, :) = 0.0d0
342  qf_p(:) = 0.0d0
343 
344  ! we suppose the same material type in the element
345  flag = gausses(1)%pMaterial%nlgeom_flag
346  elem(:, :) = ecoord(:, :)
347  totaldisp(:, :) = u(:, :)+( du(:, :)-ddu(:, :) )
348  ! Updated Lagrangian : 1st pass evaluated on the previous-iterate configuration
349  if( flag == updatelag ) elem(:, :) = ecoord(:, :)+u(:, :)+( du(:, :)-ddu(:, :) )
350  if( flag == infinitesimal ) gdispderiv(:, :) = 0.0d0
351 
352  coordsys = 0.0d0
353  coordsys(1, 1) = 1.0d0; coordsys(2, 2) = 1.0d0; coordsys(3, 3) = 1.0d0
354 
355  !========================================================
356  ! 1st pass : evaluate [ Kup ], [ Kpp ] and the pressure
357  ! residual to condense the increment ddlambda
358  !========================================================
359  do lx = 1, numofquadpoints(etype)
360 
361  call getquadpoint( etype, lx, naturalcoord(:) )
362  call getglobalderiv( etype, nn, naturalcoord, elem, det, gderiv )
363 
364  if( cdsys_id > 0 ) then
365  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
366  if( serr == -1 ) stop "Fail to setup local coordinate"
367  if( serr == -2 ) write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
368  end if
369 
370  ! Displacement gradient (B1 nonlinear term : Total Lagrangian only)
371  if( flag == totallag ) gdispderiv(1:3, 1:3) = matmul( totaldisp(1:3, 1:nn), gderiv(1:nn, 1:3) )
372 
373  if( nn_p == 1 ) spfunc_p(1) = 1.0d0
374 
375  ! Lagrange multiplier
376  lambda_lx = 0.0d0
377  do na_p = 1, nn_p
378  lambda_lx = lambda_lx+spfunc_p(na_p)*lambda(na_p)
379  end do
380 
381  mtype = gausses(lx)%pMaterial%mtype
382 
383  wg = getweight( etype, lx )*det
384 
385  ! Update strain and stress : plasticity/creep are treated as elastic in UP
386  if( iselastoplastic(gausses(lx)%pMaterial%mtype) .or. &
387  ( gausses(lx)%pMaterial%mtype == norton ) ) gausses(lx)%pMaterial%mtype = elastic
388 
389  epsth = 0.0d0
390 
391  if( present( tt ) .and. present( t0 ) ) then
392  call getshapefunc( etype, naturalcoord, spfunc )
393  ttc = dot_product( tt, spfunc )
394  tt0 = dot_product( t0, spfunc )
395  call matlmatrix_up( gausses(lx), d3, d, lambda_lx, d2, alpha_inv, g, time, tincr, coordsys, ttc )
396  if( ( iter <= 1 ) .or. ( flag == totallag ) ) then
397  ina(1) = ttc
398  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa, ierr, ina )
399  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
400  epsth(1:3) = outa(1)*( ttc-tt0 )
401  end if
402  else
403  call matlmatrix_up( gausses(lx), d3, d, lambda_lx, d2, alpha_inv, g, time, tincr, coordsys )
404  end if
405 
406  do nb = 1, nn
407  jsize1 = 3*(nb-1)+1
408  jsize2 = 3*(nb-1)+2
409  jsize3 = 3*(nb-1)+3
410  b(1, jsize1) = gderiv(nb, 1)
411  b(2, jsize1) = 0.0d0
412  b(3, jsize1) = 0.0d0
413  b(4, jsize1) = gderiv(nb, 2)
414  b(5, jsize1) = 0.0d0
415  b(6, jsize1) = gderiv(nb, 3)
416  b(1, jsize2) = 0.0d0
417  b(2, jsize2) = gderiv(nb, 2)
418  b(3, jsize2) = 0.0d0
419  b(4, jsize2) = gderiv(nb, 1)
420  b(5, jsize2) = gderiv(nb, 3)
421  b(6, jsize2) = 0.0d0
422  b(1, jsize3) = 0.0d0
423  b(2, jsize3) = 0.0d0
424  b(3, jsize3) = gderiv(nb, 3)
425  b(4, jsize3) = 0.0d0
426  b(5, jsize3) = gderiv(nb, 2)
427  b(6, jsize3) = gderiv(nb, 1)
428  end do
429 
430  ! B1 nonlinear term (Total Lagrangian only)
431  if( flag == totallag ) then
432  do nb = 1, nn
433  jsize1 = 3*(nb-1)+1
434  jsize2 = 3*(nb-1)+2
435  jsize3 = 3*(nb-1)+3
436  b1(1, jsize1) = gdispderiv(1, 1)*gderiv(nb, 1)
437  b1(2, jsize1) = gdispderiv(1, 2)*gderiv(nb, 2)
438  b1(3, jsize1) = gdispderiv(1, 3)*gderiv(nb, 3)
439  b1(4, jsize1) = gdispderiv(1, 2)*gderiv(nb, 1)+gdispderiv(1, 1)*gderiv(nb, 2)
440  b1(5, jsize1) = gdispderiv(1, 2)*gderiv(nb, 3)+gdispderiv(1, 3)*gderiv(nb, 2)
441  b1(6, jsize1) = gdispderiv(1, 3)*gderiv(nb, 1)+gdispderiv(1, 1)*gderiv(nb, 3)
442  b1(1, jsize2) = gdispderiv(2, 1)*gderiv(nb, 1)
443  b1(2, jsize2) = gdispderiv(2, 2)*gderiv(nb, 2)
444  b1(3, jsize2) = gdispderiv(2, 3)*gderiv(nb, 3)
445  b1(4, jsize2) = gdispderiv(2, 2)*gderiv(nb, 1)+gdispderiv(2, 1)*gderiv(nb, 2)
446  b1(5, jsize2) = gdispderiv(2, 2)*gderiv(nb, 3)+gdispderiv(2, 3)*gderiv(nb, 2)
447  b1(6, jsize2) = gdispderiv(2, 3)*gderiv(nb, 1)+gdispderiv(2, 1)*gderiv(nb, 3)
448  b1(1, jsize3) = gdispderiv(3, 1)*gderiv(nb, 1)
449  b1(2, jsize3) = gdispderiv(3, 2)*gderiv(nb, 2)
450  b1(3, jsize3) = gdispderiv(3, 3)*gderiv(nb, 3)
451  b1(4, jsize3) = gdispderiv(3, 2)*gderiv(nb, 1)+gdispderiv(3, 1)*gderiv(nb, 2)
452  b1(5, jsize3) = gdispderiv(3, 2)*gderiv(nb, 3)+gdispderiv(3, 3)*gderiv(nb, 2)
453  b1(6, jsize3) = gdispderiv(3, 3)*gderiv(nb, 1)+gdispderiv(3, 1)*gderiv(nb, 3)
454  end do
455  do jsize = 1, 3*nn
456  b(:, jsize) = b(:, jsize)+b1(:, jsize)
457  end do
458  end if
459 
460  do isize = 1, 3*nn
461  bd2(isize) = dot_product( b(:, isize), d2 )
462  end do
463 
464  ! [ Kup ] matrix
465  forall( isize = 1:3*nn, nb_p = 1:nn_p )
466  stiff_up(isize, nb_p) = stiff_up(isize, nb_p)+wg*bd2(isize)*spfunc_p(nb_p)
467  end forall
468 
469  ! [ Kpp ] matrix
470  forall( na_p = 1:nn_p, nb_p = 1:nn_p )
471  stiff_pp(na_p, nb_p) = stiff_pp(na_p, nb_p)-wg*alpha_inv*spfunc_p(na_p)*spfunc_p(nb_p)
472  end forall
473 
474  qf_p(1:nn_p) = qf_p(1:nn_p)+wg*spfunc_p(1:nn_p)*( g-alpha_inv*lambda_lx )
475 
476  end do ! gauss loop
477 
478  ! [ Kpp ]^-1 matrix (constant pressure: nn_p = 1)
479  if( nn_p == 1 ) then
480  stiff_pp_inv(1, 1) = 1.0d0/stiff_pp(1, 1)
481  else
482  write(6, *) 'Error: nn_p should be equal to 1.'
483  return
484  end if
485 
486  do na_p = 1, nn_p
487  stiff_up_ddu(na_p) = 0.0d0
488  do nb = 1, nn
489  do i = 1, 3
490  jsize = 3*(nb-1)+i
491  stiff_up_ddu(na_p) = stiff_up_ddu(na_p)+stiff_up(jsize, na_p)*ddu(i, nb)
492  end do
493  end do
494  end do
495 
496  do na_p = 1, nn_p
497  ddlambda(na_p) = dot_product( stiff_pp_inv(na_p, :), -qf_p-stiff_up_ddu )
498  end do
499 
500  totaldisp(:, :) = u(:, :)+du(:, :)
501  totallambda(:) = lambda(:)+ddlambda(:)
502 
503  stiff_up(:, :) = 0.0d0
504  stiff_pp(:, :) = 0.0d0
505  qf(:) = 0.0d0
506  qf_p(:) = 0.0d0
507 
508  ! Updated Lagrangian : strain increment on the mid configuration, internal
509  ! force assembled on the end (current) configuration
510  elem(:, :) = ecoord(:, :)
511  if( flag == updatelag ) then
512  elem(:, :) = ecoord(:, :)+u(:, :)+0.5d0*du(:, :)
513  elem1(:, :) = ecoord(:, :)+u(:, :)+du(:, :)
514  end if
515 
516  !========================================================
517  ! 2nd pass : update strain/stress and assemble the
518  ! condensed internal force at the new state
519  !========================================================
520  do lx = 1, numofquadpoints(etype)
521 
522  call getquadpoint( etype, lx, naturalcoord(:) )
523  call getglobalderiv( etype, nn, naturalcoord, elem, det, gderiv )
524 
525  if( cdsys_id > 0 ) then
526  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
527  if( serr == -1 ) stop "Fail to setup local coordinate"
528  if( serr == -2 ) write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
529  end if
530 
531  ! displacement gradient : step increment (UPDATELAG) or total (others)
532  if( flag == updatelag ) then
533  gdispderiv(1:3, 1:3) = matmul( du(1:3, 1:nn), gderiv(1:nn, 1:3) )
534  else
535  gdispderiv(1:3, 1:3) = matmul( totaldisp(1:3, 1:nn), gderiv(1:nn, 1:3) )
536  end if
537 
538  if( nn_p == 1 ) spfunc_p(1) = 1.0d0
539 
540  ! Lagrange multiplier
541  lambda_lx = 0.0d0
542  do na_p = 1, nn_p
543  lambda_lx = lambda_lx+spfunc_p(na_p)*totallambda(na_p)
544  end do
545 
546  mtype = gausses(lx)%pMaterial%mtype
547 
548  wg = getweight( etype, lx )*det
549 
550  ! Small (linear) strain
551  dstrain(1) = gdispderiv(1, 1)
552  dstrain(2) = gdispderiv(2, 2)
553  dstrain(3) = gdispderiv(3, 3)
554  dstrain(4) = gdispderiv(1, 2)+gdispderiv(2, 1)
555  dstrain(5) = gdispderiv(2, 3)+gdispderiv(3, 2)
556  dstrain(6) = gdispderiv(3, 1)+gdispderiv(1, 3)
557  dstrain(:) = dstrain(:)-epsth(:)
558 
559  ! Green-Lagrange strain (Total Lagrangian only)
560  if( flag == totallag ) then
561  dstrain(1) = dstrain(1)+0.5d0*dot_product( gdispderiv(:, 1), gdispderiv(:, 1) )
562  dstrain(2) = dstrain(2)+0.5d0*dot_product( gdispderiv(:, 2), gdispderiv(:, 2) )
563  dstrain(3) = dstrain(3)+0.5d0*dot_product( gdispderiv(:, 3), gdispderiv(:, 3) )
564  dstrain(4) = dstrain(4)+dot_product( gdispderiv(:, 1), gdispderiv(:, 2) )
565  dstrain(5) = dstrain(5)+dot_product( gdispderiv(:, 2), gdispderiv(:, 3) )
566  dstrain(6) = dstrain(6)+dot_product( gdispderiv(:, 1), gdispderiv(:, 3) )
567  end if
568 
569  if( flag == updatelag ) then
570  ! spin (skew part) for the objective stress rate
571  rot = 0.0d0
572  rot(1, 2) = 0.5d0*( gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
573  rot(2, 3) = 0.5d0*( gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
574  rot(1, 3) = 0.5d0*( gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
575  gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
576  else
577  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
578  end if
579 
580  if( flag == updatelag ) then
581  ! deviatoric tangent (K=0) of the increment
582  call matlmatrix_up( gausses(lx), d3, d, lambda_lx, d2, alpha_inv, g, time, tincr, coordsys )
583  ! objective (Jaumann/Hughes-Winget) update of the deviatoric stress.
584  ! pressure is imposed by the Lagrange multiplier (UP), not by the rate.
585  p_bak = ( gausses(lx)%stress_bak(1)+gausses(lx)%stress_bak(2)+gausses(lx)%stress_bak(3) )/3.0d0
586  dumstress(1, 1) = gausses(lx)%stress_bak(1)-p_bak
587  dumstress(2, 2) = gausses(lx)%stress_bak(2)-p_bak
588  dumstress(3, 3) = gausses(lx)%stress_bak(3)-p_bak
589  dumstress(1, 2) = gausses(lx)%stress_bak(4); dumstress(2, 1) = dumstress(1, 2)
590  dumstress(2, 3) = gausses(lx)%stress_bak(5); dumstress(3, 2) = dumstress(2, 3)
591  dumstress(3, 1) = gausses(lx)%stress_bak(6); dumstress(1, 3) = dumstress(3, 1)
592  trd = dstrain(1)+dstrain(2)+dstrain(3)
593  dum(:, :) = dumstress+matmul( rot, dumstress )-matmul( dumstress, rot )-dumstress*trd
594  dstress(1:6) = matmul( d(1:6, 1:6), dstrain(1:6) )
595  gausses(lx)%stress(1) = dum(1, 1)+dstress(1)+lambda_lx
596  gausses(lx)%stress(2) = dum(2, 2)+dstress(2)+lambda_lx
597  gausses(lx)%stress(3) = dum(3, 3)+dstress(3)+lambda_lx
598  gausses(lx)%stress(4) = dum(1, 2)+dstress(4)
599  gausses(lx)%stress(5) = dum(2, 3)+dstress(5)
600  gausses(lx)%stress(6) = dum(3, 1)+dstress(6)
601  else
602  call stressupdate_up( gausses(lx), d3, dstrain, gausses(lx)%stress, lambda_lx, g, coordsys )
603  end if
604 
605  gausses(lx)%stress_out(1:6) = gausses(lx)%stress(1:6)
606  gausses(lx)%strain_out(1:6) = gausses(lx)%strain(1:6)
607 
608  epsth = 0.0d0
609 
610  if( present( tt ) .and. present( t0 ) ) then
611  call getshapefunc( etype, naturalcoord, spfunc )
612  ttc = dot_product( tt, spfunc )
613  tt0 = dot_product( t0, spfunc )
614  call matlmatrix_up( gausses(lx), d3, d, lambda_lx, d2, alpha_inv, g, time, tincr, coordsys, ttc )
615  if( ( iter <= 1 ) .or. ( flag == totallag ) ) then
616  ina(1) = ttc
617  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa, ierr, ina )
618  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
619  epsth(1:3) = outa(1)*( ttc-tt0 )
620  end if
621  else
622  call matlmatrix_up( gausses(lx), d3, d, lambda_lx, d2, alpha_inv, g, time, tincr, coordsys )
623  end if
624 
625  ! Updated Lagrangian : switch to the end (current) configuration for the
626  ! B-matrix and the internal force / coupling assembly
627  if( flag == updatelag ) then
628  call getglobalderiv( etype, nn, naturalcoord, elem1, det, gderiv )
629  wg = getweight( etype, lx )*det
630  end if
631 
632  do nb = 1, nn
633  jsize1 = 3*(nb-1)+1
634  jsize2 = 3*(nb-1)+2
635  jsize3 = 3*(nb-1)+3
636  b(1, jsize1) = gderiv(nb, 1)
637  b(2, jsize1) = 0.0d0
638  b(3, jsize1) = 0.0d0
639  b(4, jsize1) = gderiv(nb, 2)
640  b(5, jsize1) = 0.0d0
641  b(6, jsize1) = gderiv(nb, 3)
642  b(1, jsize2) = 0.0d0
643  b(2, jsize2) = gderiv(nb, 2)
644  b(3, jsize2) = 0.0d0
645  b(4, jsize2) = gderiv(nb, 1)
646  b(5, jsize2) = gderiv(nb, 3)
647  b(6, jsize2) = 0.0d0
648  b(1, jsize3) = 0.0d0
649  b(2, jsize3) = 0.0d0
650  b(3, jsize3) = gderiv(nb, 3)
651  b(4, jsize3) = 0.0d0
652  b(5, jsize3) = gderiv(nb, 2)
653  b(6, jsize3) = gderiv(nb, 1)
654  end do
655 
656  ! B1 nonlinear term (Total Lagrangian only)
657  if( flag == totallag ) then
658  do nb = 1, nn
659  jsize1 = 3*(nb-1)+1
660  jsize2 = 3*(nb-1)+2
661  jsize3 = 3*(nb-1)+3
662  b1(1, jsize1) = gdispderiv(1, 1)*gderiv(nb, 1)
663  b1(2, jsize1) = gdispderiv(1, 2)*gderiv(nb, 2)
664  b1(3, jsize1) = gdispderiv(1, 3)*gderiv(nb, 3)
665  b1(4, jsize1) = gdispderiv(1, 2)*gderiv(nb, 1)+gdispderiv(1, 1)*gderiv(nb, 2)
666  b1(5, jsize1) = gdispderiv(1, 2)*gderiv(nb, 3)+gdispderiv(1, 3)*gderiv(nb, 2)
667  b1(6, jsize1) = gdispderiv(1, 3)*gderiv(nb, 1)+gdispderiv(1, 1)*gderiv(nb, 3)
668  b1(1, jsize2) = gdispderiv(2, 1)*gderiv(nb, 1)
669  b1(2, jsize2) = gdispderiv(2, 2)*gderiv(nb, 2)
670  b1(3, jsize2) = gdispderiv(2, 3)*gderiv(nb, 3)
671  b1(4, jsize2) = gdispderiv(2, 2)*gderiv(nb, 1)+gdispderiv(2, 1)*gderiv(nb, 2)
672  b1(5, jsize2) = gdispderiv(2, 2)*gderiv(nb, 3)+gdispderiv(2, 3)*gderiv(nb, 2)
673  b1(6, jsize2) = gdispderiv(2, 3)*gderiv(nb, 1)+gdispderiv(2, 1)*gderiv(nb, 3)
674  b1(1, jsize3) = gdispderiv(3, 1)*gderiv(nb, 1)
675  b1(2, jsize3) = gdispderiv(3, 2)*gderiv(nb, 2)
676  b1(3, jsize3) = gdispderiv(3, 3)*gderiv(nb, 3)
677  b1(4, jsize3) = gdispderiv(3, 2)*gderiv(nb, 1)+gdispderiv(3, 1)*gderiv(nb, 2)
678  b1(5, jsize3) = gdispderiv(3, 2)*gderiv(nb, 3)+gdispderiv(3, 3)*gderiv(nb, 2)
679  b1(6, jsize3) = gdispderiv(3, 3)*gderiv(nb, 1)+gdispderiv(3, 1)*gderiv(nb, 3)
680  end do
681  do jsize = 1, 3*nn
682  b(:, jsize) = b(:, jsize)+b1(:, jsize)
683  end do
684  end if
685 
686  do isize = 1, 3*nn
687  bd2(isize) = dot_product( b(:, isize), d2 )
688  end do
689 
690  ! [ Kup ] matrix
691  forall( isize = 1:3*nn, nb_p = 1:nn_p )
692  stiff_up(isize, nb_p) = stiff_up(isize, nb_p)+wg*bd2(isize)*spfunc_p(nb_p)
693  end forall
694 
695  ! [ Kpp ] matrix
696  forall( na_p = 1:nn_p, nb_p = 1:nn_p )
697  stiff_pp(na_p, nb_p) = stiff_pp(na_p, nb_p)-wg*alpha_inv*spfunc_p(na_p)*spfunc_p(nb_p)
698  end forall
699 
700  ! calculate the Internal Force
701  qf(1:3*nn) = qf(1:3*nn)+wg*matmul( gausses(lx)%stress(1:6), b(1:6, 1:3*nn) )
702 
703  qf_p(1:nn_p) = qf_p(1:nn_p)+wg*spfunc_p(1:nn_p)*( g-alpha_inv*lambda_lx )
704 
705  end do ! gauss loop
706 
707  ! [ Kpp ]^-1 matrix (constant pressure: nn_p = 1)
708  if( nn_p == 1 ) then
709  stiff_pp_inv(1, 1) = 1.0d0/stiff_pp(1, 1)
710  else
711  write(6, *) 'Error: nn_p should be equal to 1.'
712  return
713  end if
714 
715  stiff_up_stiff_pp_inv(1:3*nn, 1:nn_p) = matmul( stiff_up(1:3*nn, 1:nn_p), stiff_pp_inv(1:nn_p, 1:nn_p) )
716 
717  do isize = 1, 3*nn
718  stiff_up_stiff_pp_inv_qf_p(isize) = dot_product( stiff_up_stiff_pp_inv(isize, :), qf_p )
719  end do
720 
721  ! condensed internal force : qf - [ Kup ] [ Kpp ]^-1 qf_p
722  do isize = 1, 3*nn
723  qf(isize) = qf(isize)-stiff_up_stiff_pp_inv_qf_p(isize)
724  end do
725 
726  end subroutine update_c3_up
727 
728 end module m_static_lib_3d_up
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::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_matmatrix
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
m_matmatrix::stressupdate_up
subroutine stressupdate_up(gauss, sectType, strain, stress, lambda, g, cdsys)
Update stress for the u-p mixed formulation. Deviatoric stress from the material law (K=0) plus the p...
Definition: calMatMatrix.f90:450
m_matmatrix::matlmatrix_up
subroutine matlmatrix_up(gauss, sectType, D, lambda, d2, alpha_inv, g, time, dtime, cdsys, temperature)
Constitutive matrix for the u-p mixed formulation. Returns the deviatoric tangent D (K=0),...
Definition: calMatMatrix.f90:411
m_static_lib_3d_up
This module provides functions of u-p mixed (UP) solid elements.
Definition: static_LIB_3d_up.f90:6
m_static_lib_3d_up::update_c3_up
subroutine update_c3_up(etype, nn, ecoord, u, du, ddu, cdsys_ID, coords, qf, gausses, iter, time, tincr, nn_p, lambda, ddlambda, tt, t0)
Update strain and stress inside u-p mixed (UP) solid element.
Definition: static_LIB_3d_up.f90:282
m_static_lib_3d_up::stf_c3_up
subroutine stf_c3_up(etype, nn, ecoord, gausses, stiff, cdsys_ID, coords, time, tincr, nn_p, lambda, u, temperature)
This subroutine calculate stiff matrix of u-p mixed solid elements.
Definition: static_LIB_3d_up.f90:19
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_utilities
This module provides aux functions.
Definition: utilities.f90:6
mhyperelastic
This module provides functions for hyperelastic calculation.
Definition: Hyperelastic.f90:6
mmaterial
This module summarizes all information of material properties.
Definition: material.f90:6
mmaterial::m_exapnsion
integer(kind=kint), parameter m_exapnsion
Definition: material.f90:105
mmaterial::mc_themoexp
character(len=dict_key_length) mc_themoexp
Definition: material.f90:143
mmaterial::d3
integer(kind=kint), parameter d3
Definition: material.f90:84
mmaterial::elastic
integer(kind=kint), parameter elastic
Definition: material.f90:65
mmaterial::totallag
integer(kind=kint), parameter totallag
Definition: material.f90:14
mmaterial::norton
integer(kind=kint), parameter norton
Definition: material.f90:78
mmaterial::infinitesimal
integer(kind=kint), parameter infinitesimal
Definition: material.f90:13
mmaterial::updatelag
integer(kind=kint), parameter updatelag
Definition: material.f90:15
mmaterial::iselastoplastic
logical function iselastoplastic(mtype)
If it is an elastoplastic material?
Definition: material.f90:363
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:16