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