FrontISTR  5.7.1
Large-scale structural analysis program with finit element method
static_LIB_3d.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 !-------------------------------------------------------------------------------
7 
8  use hecmw, only : kint, kreal
9  use elementinfo
10 
11  implicit none
12 
13  real(kind=kreal), parameter, private :: i3(3,3) = reshape((/ &
14  & 1.d0, 0.d0, 0.d0, 0.d0, 1.d0, 0.d0, 0.d0, 0.d0, 1.d0/), (/3,3/))
15 
16 contains
17  !----------------------------------------------------------------------*
18  subroutine geomat_c3(stress, mat)
19  !----------------------------------------------------------------------*
20 
21  real(kind=kreal), intent(in) :: stress(6)
22  real(kind=kreal), intent(out) :: mat(6, 6)
23 
24  !---------------------------------------------------------------------
25 
26  mat(1, 1) = 2.0d0*stress(1); mat(1, 2) = 0.0d0; mat(1, 3) = 0.0d0
27  mat(1, 4) = stress(4); mat(1, 5) = 0.0d0; mat(1, 6) = stress(6)
28  mat(2, 1) = mat(1, 2); mat(2, 2) = 2.0d0*stress(2); mat(2, 3) = 0.0d0
29  mat(2, 4) = stress(4); mat(2, 5) = stress(5); mat(2, 6) = 0.0d0
30  mat(3, 1) = mat(1, 3); mat(3, 2) = mat(2, 3); mat(3, 3) = 2.0d0*stress(3)
31  mat(3, 4) = 0.0d0; mat(3, 5) = stress(5); mat(3, 6) = stress(6)
32 
33  mat(4, 1) = mat(1, 4); mat(4, 2) = mat(2, 4); mat(4, 3) = mat(3, 4)
34  mat(4, 4) = 0.5d0*( stress(1)+stress(2) ); mat(4, 5) = 0.5d0*stress(6); mat(4, 6) = 0.5d0*stress(5)
35  mat(5, 1) = mat(1, 5); mat(5, 2) = mat(2, 5); mat(5, 3) = mat(3, 5)
36  mat(5, 4) = mat(4, 5); mat(5, 5) = 0.5d0*( stress(3)+stress(2) ); mat(5, 6) = 0.5d0*stress(4)
37  mat(6, 1) = mat(1, 6); mat(6, 2) = mat(2, 6); mat(6, 3) = mat(3, 6)
38  mat(6, 4) = mat(4, 6); mat(6, 5) = mat(5, 6); mat(6, 6) = 0.5d0*( stress(1)+stress(3) );
39 
40  end subroutine
41 
42 
43  !=====================================================================*
45  !
49  !----------------------------------------------------------------------*
50  subroutine stf_c3 &
51  (etype, nn, ecoord, gausses, stiff, cdsys_id, coords, &
52  time, tincr, u ,temperature)
53  !----------------------------------------------------------------------*
54 
55  use mmechgauss
56  use m_matmatrix
57  use m_common_struct
58 
59  !---------------------------------------------------------------------
60 
61  integer(kind=kint), intent(in) :: etype
62  integer(kind=kint), intent(in) :: nn
63  real(kind=kreal), intent(in) :: ecoord(3,nn)
64  type(tgaussstatus), intent(in) :: gausses(:)
65  real(kind=kreal), intent(out) :: stiff(:,:)
66  integer(kind=kint), intent(in) :: cdsys_id
67  real(kind=kreal), intent(inout) :: coords(3,3)
68  real(kind=kreal), intent(in) :: time
69  real(kind=kreal), intent(in) :: tincr
70  real(kind=kreal), intent(in) :: temperature(nn)
71  real(kind=kreal), intent(in), optional :: u(:,:)
72 
73  !---------------------------------------------------------------------
74 
75  integer(kind=kint) :: flag
76  integer(kind=kint), parameter :: ndof = 3
77  real(kind=kreal) :: d(6, 6), b(6, ndof*nn), db(6, ndof*nn)
78  real(kind=kreal) :: gderiv(nn, 3), stress(6), mat(6, 6)
79  real(kind=kreal) :: det, wg
80  integer(kind=kint) :: i, j, lx, serr
81  real(kind=kreal) :: temp, naturalcoord(3)
82  real(kind=kreal) :: spfunc(nn), gdispderiv(3, 3)
83  real(kind=kreal) :: b1(6, ndof*nn), coordsys(3, 3)
84  real(kind=kreal) :: smat(9, 9), elem(3, nn)
85  real(kind=kreal) :: bn(9, ndof*nn), sbn(9, ndof*nn)
86 
87  !---------------------------------------------------------------------
88 
89  stiff(:, :) = 0.0d0
90  ! we suppose the same material type in the element
91  flag = gausses(1)%pMaterial%nlgeom_flag
92  if( .not. present(u) ) flag = infinitesimal ! enforce to infinitesimal deformation analysis
93  elem(:, :) = ecoord(:, :)
94  if( flag == updatelag ) elem(:, :) = ecoord(:, :)+u(:, :)
95 
96  do lx = 1, numofquadpoints(etype)
97 
98  call getquadpoint( etype, lx, naturalcoord(:) )
99  call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
100 
101  if( cdsys_id > 0 ) then
102  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:, :), serr )
103  if( serr == -1 ) stop "Fail to setup local coordinate"
104  if( serr == -2 ) then
105  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
106  end if
107  end if
108 
109  call getshapefunc(etype, naturalcoord, spfunc)
110  temp = dot_product(temperature, spfunc)
111  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
112 
113  if( flag == updatelag ) then
114  call geomat_c3( gausses(lx)%stress, mat )
115  d(:, :) = d(:, :)-mat
116  endif
117 
118  wg = getweight(etype, lx)*det
119  b(1:6, 1:nn*ndof) = 0.0d0
120  do j = 1, nn
121  b(1, 3*j-2)=gderiv(j, 1)
122  b(2, 3*j-1)=gderiv(j, 2)
123  b(3, 3*j )=gderiv(j, 3)
124  b(4, 3*j-2)=gderiv(j, 2)
125  b(4, 3*j-1)=gderiv(j, 1)
126  b(5, 3*j-1)=gderiv(j, 3)
127  b(5, 3*j )=gderiv(j, 2)
128  b(6, 3*j-2)=gderiv(j, 3)
129  b(6, 3*j )=gderiv(j, 1)
130  enddo
131 
132  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
133  if( flag == totallag ) then
134  ! ---dudx(i, j) ==> gdispderiv(i, j)
135  gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
136  b1(1:6, 1:nn*ndof)=0.0d0
137  do j=1, nn
138  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
139  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
140  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
141  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
142  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
143  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
144  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
145  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
146  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
147  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
148  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
149  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
150  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
151  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
152  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
153  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
154  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
155  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
156  end do
157  ! ---BL = BL0 + BL1
158  do j=1, nn*ndof
159  b(:, j) = b(:, j)+b1(:, j)
160  end do
161  end if
162 
163  db(1:6, 1:nn*ndof) = matmul( d, b(1:6, 1:nn*ndof) )
164  do j=1,nn*ndof
165  do i=1,nn*ndof
166  stiff(i, j) = stiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
167  enddo
168  enddo
169 
170  ! calculate the stress matrix ( TOTAL LAGRANGE METHOD )
171  if( flag == totallag .OR. flag==updatelag ) then
172  stress(1:6) = gausses(lx)%stress
173  bn(1:9, 1:nn*ndof) = 0.0d0
174  do j = 1, nn
175  bn(1, 3*j-2) = gderiv(j, 1)
176  bn(2, 3*j-1) = gderiv(j, 1)
177  bn(3, 3*j ) = gderiv(j, 1)
178  bn(4, 3*j-2) = gderiv(j, 2)
179  bn(5, 3*j-1) = gderiv(j, 2)
180  bn(6, 3*j ) = gderiv(j, 2)
181  bn(7, 3*j-2) = gderiv(j, 3)
182  bn(8, 3*j-1) = gderiv(j, 3)
183  bn(9, 3*j ) = gderiv(j, 3)
184  end do
185  smat(:, :) = 0.0d0
186  do j = 1, 3
187  smat(j , j ) = stress(1)
188  smat(j , j+3) = stress(4)
189  smat(j , j+6) = stress(6)
190  smat(j+3, j ) = stress(4)
191  smat(j+3, j+3) = stress(2)
192  smat(j+3, j+6) = stress(5)
193  smat(j+6, j ) = stress(6)
194  smat(j+6, j+3) = stress(5)
195  smat(j+6, j+6) = stress(3)
196  end do
197  sbn(1:9, 1:nn*ndof) = matmul( smat(1:9, 1:9), bn(1:9, 1:nn*ndof) )
198  do j=1,nn*ndof
199  do i=1,nn*ndof
200  stiff(i, j) = stiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
201  enddo
202  enddo
203 
204  end if
205 
206  enddo ! gauss roop
207 
208  end subroutine stf_c3
209 
210 
212  !----------------------------------------------------------------------*
213  subroutine dl_c3(etype, nn, XX, YY, ZZ, RHO, LTYPE, PARAMS, VECT, NSIZE)
214  !--------------------------------------------------------------------*
215  !**
216  !** SET DLOAD
217  !**
218  ! BX LTYPE=1 :BODY FORCE IN X-DIRECTION
219  ! BY LTYPE=2 :BODY FORCE IN Y-DIRECTION
220  ! BZ LTYPE=3 :BODY FORCE IN Z-DIRECTION
221  ! GRAV LTYPE=4 :GRAVITY FORCE
222  ! CENT LTYPE=5 :CENTRIFUGAL LOAD
223  ! P1 LTYPE=10 :TRACTION IN NORMAL-DIRECTION FOR FACE-1
224  ! P2 LTYPE=20 :TRACTION IN NORMAL-DIRECTION FOR FACE-2
225  ! P3 LTYPE=30 :TRACTION IN NORMAL-DIRECTION FOR FACE-3
226  ! P4 LTYPE=40 :TRACTION IN NORMAL-DIRECTION FOR FACE-4
227  ! P5 LTYPE=50 :TRACTION IN NORMAL-DIRECTION FOR FACE-5
228  ! P6 LTYPE=60 :TRACTION IN NORMAL-DIRECTION FOR FACE-6
229  ! I/F VARIABLES
230  integer(kind=kint), intent(in) :: etype, nn
231  real(kind=kreal), intent(in) :: xx(:),yy(:),zz(:)
232  real(kind=kreal), intent(in) :: params(0:6)
233  real(kind=kreal), intent(inout) :: vect(:)
234  real(kind=kreal) rho
235  integer(kind=kint) ltype,nsize
236 
237  ! LOCAL VARIABLES
238  integer(kind=kint) ndof
239  parameter(ndof=3)
240  real(kind=kreal) h(nn)
241  real(kind=kreal) plx(nn), ply(nn), plz(nn)
242  real(kind=kreal) xj(3, 3), det, wg
243  integer(kind=kint) ivol, isuf
244  integer(kind=kint) nod(nn)
245  integer(kind=kint) IG2, LX, I , SURTYPE, NSUR
246  real(kind=kreal) vx, vy, vz, xcod, ycod, zcod
247  real(kind=kreal) ax, ay, az, rx, ry, rz, hx, hy, hz, val
248  real(kind=kreal) phx, phy, phz
249  real(kind=kreal) coefx, coefy, coefz
250  real(kind=kreal) normal(3), localcoord(3), elecoord(3, nn), deriv(nn, 3)
251 
252  ax = 0.0d0; ay = 0.0d0; az = 0.0d0; rx = 0.0d0; ry = 0.0d0; rz = 0.0d0;
253  !
254  ! SET VALUE
255  !
256  val = params(0)
257  !
258  ! SELECTION OF LOAD TYPE
259  !
260  ivol=0
261  isuf=0
262  if( ltype.LT.10 ) then
263  ivol=1
264  if( ltype.EQ.5 ) then
265  ax=params(1)
266  ay=params(2)
267  az=params(3)
268  rx=params(4)
269  ry=params(5)
270  rz=params(6)
271  endif
272  else if( ltype.GE.10 ) then
273  isuf=1
274  call getsubface( etype, ltype/10, surtype, nod )
275  nsur = getnumberofnodes( surtype )
276  endif
277  ! CLEAR VECT
278  nsize=nn*ndof
279  vect(1:nsize)=0.0d0
280  !** SURFACE LOAD
281  if( isuf==1 ) then
282  ! INTEGRATION OVER SURFACE
283  do i=1,nsur
284  elecoord(1,i)=xx(nod(i))
285  elecoord(2,i)=yy(nod(i))
286  elecoord(3,i)=zz(nod(i))
287  enddo
288  do ig2=1,numofquadpoints( surtype )
289  call getquadpoint( surtype, ig2, localcoord(1:2) )
290  call getshapefunc( surtype, localcoord(1:2), h(1:nsur) )
291 
292  wg=getweight( surtype, ig2 )
293  normal=surfacenormal( surtype, nsur, localcoord(1:2), elecoord(:,1:nsur) )
294  do i=1,nsur
295  vect(3*nod(i)-2)=vect(3*nod(i)-2)+val*wg*h(i)*normal(1)
296  vect(3*nod(i)-1)=vect(3*nod(i)-1)+val*wg*h(i)*normal(2)
297  vect(3*nod(i) )=vect(3*nod(i) )+val*wg*h(i)*normal(3)
298  enddo
299  enddo
300  endif
301  !** VOLUME LOAD
302  if( ivol==1 ) then
303  plx(:)=0.0d0
304  ply(:)=0.0d0
305  plz(:)=0.0d0
306  ! LOOP FOR INTEGRATION POINTS
307  do lx=1,numofquadpoints( etype )
308  call getquadpoint( etype, lx, localcoord )
309  call getshapefunc( etype, localcoord, h(1:nn) )
310  call getshapederiv( etype, localcoord, deriv )
311  ! JACOBI MATRIX
312  xj(1,1:3)= matmul( xx(1:nn), deriv(1:nn,1:3) )
313  xj(2,1:3)= matmul( yy(1:nn), deriv(1:nn,1:3) )
314  xj(3,1:3)= matmul( zz(1:nn), deriv(1:nn,1:3) )
315  !DETERMINANT OF JACOBIAN
316  det=xj(1,1)*xj(2,2)*xj(3,3) &
317  +xj(2,1)*xj(3,2)*xj(1,3) &
318  +xj(3,1)*xj(1,2)*xj(2,3) &
319  -xj(3,1)*xj(2,2)*xj(1,3) &
320  -xj(2,1)*xj(1,2)*xj(3,3) &
321  -xj(1,1)*xj(3,2)*xj(2,3)
322 
323  coefx=1.0
324  coefy=1.0
325  coefz=1.0
326  ! CENTRIFUGAL LOAD
327  if( ltype==5 ) then
328  xcod=dot_product( h(1:nn),xx(1:nn) )
329  ycod=dot_product( h(1:nn),yy(1:nn) )
330  zcod=dot_product( h(1:nn),zz(1:nn) )
331  hx=ax+((xcod-ax)*rx+(ycod-ay)*ry+(zcod-az)*rz)/(rx**2+ry**2+rz**2)*rx
332  hy=ay+((xcod-ax)*rx+(ycod-ay)*ry+(zcod-az)*rz)/(rx**2+ry**2+rz**2)*ry
333  hz=az+((xcod-ax)*rx+(ycod-ay)*ry+(zcod-az)*rz)/(rx**2+ry**2+rz**2)*rz
334  phx=xcod-hx
335  phy=ycod-hy
336  phz=zcod-hz
337  coefx=rho*val*dabs(val)*phx
338  coefy=rho*val*dabs(val)*phy
339  coefz=rho*val*dabs(val)*phz
340  end if
341 
342  wg=getweight( etype, lx )*det
343  do i=1,nn
344  plx(i)=plx(i)+h(i)*wg*coefx
345  ply(i)=ply(i)+h(i)*wg*coefy
346  plz(i)=plz(i)+h(i)*wg*coefz
347  enddo
348  enddo
349  if( ltype.EQ.1) then
350  do i=1,nn
351  vect(3*i-2)=val*plx(i)
352  enddo
353  else if( ltype.EQ.2 ) then
354  do i=1,nn
355  vect(3*i-1)=val*ply(i)
356  enddo
357  else if( ltype.EQ.3 ) then
358  do i=1,nn
359  vect(3*i )=val*plz(i)
360  enddo
361  else if( ltype.EQ.4 ) then
362  vx=params(1)
363  vy=params(2)
364  vz=params(3)
365  vx=vx/sqrt(params(1)**2+params(2)**2+params(3)**2)
366  vy=vy/sqrt(params(1)**2+params(2)**2+params(3)**2)
367  vz=vz/sqrt(params(1)**2+params(2)**2+params(3)**2)
368  do i=1,nn
369  vect(3*i-2)=val*plx(i)*rho*vx
370  vect(3*i-1)=val*ply(i)*rho*vy
371  vect(3*i )=val*plz(i)*rho*vz
372  enddo
373  else if( ltype.EQ.5 ) then
374  do i=1,nn
375  vect(3*i-2)=plx(i)
376  vect(3*i-1)=ply(i)
377  vect(3*i )=plz(i)
378  enddo
379  end if
380  endif
381 
382  end subroutine dl_c3
383 
385  !----------------------------------------------------------------------*
386  subroutine tload_c3 &
387  (etype, nn, xx, yy, zz, tt, t0, gausses, &
388  vect, cdsys_id, coords)
389  !----------------------------------------------------------------------*
390 
391  use m_fstr
392  use mmechgauss
393  use m_matmatrix
394  use m_utilities
395 
396  !---------------------------------------------------------------------
397 
398  integer(kind=kint), parameter :: ndof = 3
399  integer(kind=kint), intent(in) :: etype, nn
400  type(tgaussstatus), intent(in) :: gausses(:)
401  real(kind=kreal), intent(in) :: xx(nn), yy(nn), zz(nn)
402  real(kind=kreal), intent(in) :: tt(nn),t0(nn)
403  real(kind=kreal), intent(out) :: vect(nn*ndof)
404  integer(kind=kint), intent(in) :: cdsys_ID
405  real(kind=kreal), intent(inout) :: coords(3, 3)
406 
407  !---------------------------------------------------------------------
408 
409  real(kind=kreal) :: alp, alp0, d(6, 6), b(6, ndof*nn)
410  real(kind=kreal) :: det, ecoord(3, nn)
411  integer(kind=kint) :: j, LX, serr
412  real(kind=kreal) :: epsth(6),sgm(6), h(nn), alpo(3), alpo0(3), coordsys(3, 3)
413  real(kind=kreal) :: naturalcoord(3), gderiv(nn, 3)
414  real(kind=kreal) :: wg, outa(1), ina(1)
415  real(kind=kreal) :: tempc, temp0, thermal_eps, tm(6, 6)
416  logical :: ierr, matlaniso
417 
418  !---------------------------------------------------------------------
419 
420  matlaniso = .false.
421 
422  if( cdsys_id > 0 ) then ! cannot define aniso expansion when no local coord defined
423  ina = tt(1)
424  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
425  if( .not. ierr ) matlaniso = .true.
426  end if
427 
428  vect(:) = 0.0d0
429 
430  ecoord(1, :) = xx(:)
431  ecoord(2, :) = yy(:)
432  ecoord(3, :) = zz(:)
433 
434  ! LOOP FOR INTEGRATION POINTS
435  do lx = 1, numofquadpoints(etype)
436 
437  call getquadpoint( etype, lx, naturalcoord(:) )
438  call getshapefunc( etype, naturalcoord, h(1:nn) )
439  call getglobalderiv( etype, nn, naturalcoord, ecoord, det, gderiv )
440 
441  if( cdsys_id > 0 ) then
442  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys, serr )
443  if( serr == -1 ) stop "Fail to setup local coordinate"
444  if( serr == -2 ) then
445  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
446  end if
447  end if
448 
449  ! WEIGHT VALUE AT GAUSSIAN POINT
450  wg = getweight(etype, lx)*det
451  b(1:6,1:nn*ndof)=0.0d0
452  do j=1,nn
453  b(1,3*j-2)=gderiv(j,1)
454  b(2,3*j-1)=gderiv(j,2)
455  b(3,3*j )=gderiv(j,3)
456  b(4,3*j-2)=gderiv(j,2)
457  b(4,3*j-1)=gderiv(j,1)
458  b(5,3*j-1)=gderiv(j,3)
459  b(5,3*j )=gderiv(j,2)
460  b(6,3*j-2)=gderiv(j,3)
461  b(6,3*j )=gderiv(j,1)
462  enddo
463 
464  tempc = dot_product( h(1:nn), tt(1:nn) )
465  temp0 = dot_product( h(1:nn), t0(1:nn) )
466 
467  call matlmatrix( gausses(lx), d3, d, 1.d0, 0.0d0, coordsys, tempc )
468 
469  ina(1) = tempc
470  if( matlaniso ) then
471  call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo(:), ierr, ina )
472  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
473  else
474  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
475  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
476  alp = outa(1)
477  end if
478  ina(1) = temp0
479  if( matlaniso ) then
480  call fetch_tabledata( mc_orthoexp, gausses(lx)%pMaterial%dict, alpo0(:), ierr, ina )
481  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
482  else
483  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
484  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
485  alp0 = outa(1)
486  end if
487 
488  !**
489  !** THERMAL strain
490  !**
491  if( matlaniso ) then
492  do j=1,3
493  epsth(j) = alpo(j)*(tempc-ref_temp)-alpo0(j)*(temp0-ref_temp)
494  end do
495  epsth(4:6) = 0.0d0
496  call transformation(coordsys, tm)
497  epsth(:) = matmul( epsth(:), tm ) ! to global coord
498  epsth(4:6) = epsth(4:6)*2.0d0
499  else
500  thermal_eps=alp*(tempc-ref_temp)-alp0*(temp0-ref_temp)
501  epsth(1:3) = thermal_eps
502  epsth(4:6) = 0.0d0
503  end if
504 
505  !**
506  !** SET SGM {s}=[D]{e}
507  !**
508  sgm(:) = matmul( d(:, :), epsth(:) )
509  !**
510  !** CALCULATE LOAD {F}=[B]T{e}
511  !**
512  vect(1:nn*ndof) = vect(1:nn*ndof)+matmul( sgm(:),b(:, :) )*wg
513 
514  end do
515 
516  end subroutine tload_c3
517 
518  subroutine cal_thermal_expansion_c3( tt0, ttc, material, coordsys, matlaniso, EPSTH )
519  use m_fstr
520  use m_utilities
521  real(kind=kreal), INTENT(IN) :: tt0
522  real(kind=kreal), INTENT(IN) :: ttc
523  type( tmaterial ), INTENT(IN) :: material
524  real(kind=kreal), INTENT(IN) :: coordsys(3,3)
525  logical, INTENT(IN) :: matlaniso
526  real(kind=kreal), INTENT(OUT) :: epsth(6)
527 
528  integer(kind=kint) :: j
529  real(kind=kreal) :: ina(1), outa(1)
530  logical :: ierr
531  real(kind=kreal) :: alp, alp0, alpo(3), alpo0(3), tm(6,6)
532 
533  epsth = 0.d0
534  if( dabs(ttc-tt0) < 1.d-14 ) return
535 
536  ina(1) = ttc
537  if( matlaniso ) then
538  call fetch_tabledata( mc_orthoexp, material%dict, alpo(:), ierr, ina )
539  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
540  else
541  call fetch_tabledata( mc_themoexp, material%dict, outa(:), ierr, ina )
542  if( ierr ) outa(1) = material%variables(m_exapnsion)
543  alp = outa(1)
544  end if
545  ina(1) = tt0
546  if( matlaniso ) then
547  call fetch_tabledata( mc_orthoexp, material%dict, alpo0(:), ierr, ina )
548  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
549  else
550  call fetch_tabledata( mc_themoexp, material%dict, outa(:), ierr, ina )
551  if( ierr ) outa(1) = material%variables(m_exapnsion)
552  alp0 = outa(1)
553  end if
554  if( matlaniso ) then
555  do j=1,3
556  epsth(j) = alpo(j)*(ttc-ref_temp)-alpo0(j)*(tt0-ref_temp)
557  end do
558  call transformation( coordsys(:, :), tm)
559  epsth(:) = matmul( epsth(:), tm ) ! to global coord
560  epsth(4:6) = epsth(4:6)*2.0d0
561  else
562  epsth(1:3)=alp*(ttc-ref_temp)-alp0*(tt0-ref_temp)
563  end if
564  end subroutine
565 
566  !Hughes, T. J. R., and J. Winget, ?gFinite Rotation Effects in Numerical Integration of
567  ! Rate Constitutive Equations Arising in Large Deformation Analysis,?h
568  ! International Journal for Numerical Methods in Engineering, vol. 15, pp. 1862-1867, 1980.
569  subroutine hughes_winget_rotation_3d( rot, stress_in, stress_out )
570  use m_utilities
571  real(kind=kreal), intent(in) :: rot(3,3)
572  real(kind=kreal), intent(in) :: stress_in(3,3)
573  real(kind=kreal), intent(out) :: stress_out(3,3)
574 
575  real(kind=kreal) :: dr(3,3)
576 
577  !calc dR=(I-0.5*dW)^-1*(I+0.5*dW)
578  dr = i3-0.5d0*rot
579  call calinverse(3, dr)
580  dr = matmul(dr,i3+0.5d0*rot)
581 
582  stress_out = matmul(dr,stress_in)
583  stress_out = matmul(stress_out,transpose(dr))
584  end subroutine
585 
586  subroutine update_stress3d( flag, gauss, rot, dstrain, F, coordsys, time, tincr, ttc, tt0, ttn, hdflag )
587  use m_fstr
588  use m_matmatrix
589  use mmechgauss
590  use m_utilities
591 
592  type(tgaussstatus), intent(inout) :: gauss
593  integer(kind=kint), intent(in) :: flag
594  real(kind=kreal), intent(in) :: rot(3,3)
595  real(kind=kreal), intent(in) :: dstrain(6)
596  real(kind=kreal), intent(in) :: f(3,3) !deformation gradient (used for ss_out)
597  real(kind=kreal), intent(in) :: coordsys(3,3)
598  real(kind=kreal), intent(in) :: time
599  real(kind=kreal), intent(in) :: tincr
600  real(kind=kreal), intent(in) :: ttc
601  real(kind=kreal), intent(in) :: tt0
602  real(kind=kreal), intent(in) :: ttn
603  integer(kind=kint), intent(in), optional :: hdflag
604 
605  integer(kind=kint) :: mtype, i, j, k
606  integer(kind=kint) :: isep, hdflag_in
607  real(kind=kreal) :: d(6,6), dstress(6), dumstress(3,3), dum(3,3), trd, det
608  real(kind=kreal) :: tensor(6)
609  real(kind=kreal) :: eigval(3)
610  real(kind=kreal) :: princ(3,3), norm
611 
612  mtype = gauss%pMaterial%mtype
613 
614  if( iselastoplastic(mtype) .OR. mtype == norton )then
615  isep = 1
616  else
617  isep = 0
618  endif
619 
620  hdflag_in = 0
621  if( present(hdflag) ) hdflag_in = hdflag
622 
623  call matlmatrix( gauss, d3, d, time, tincr, coordsys, ttc, isep, hdflag=hdflag_in )
624 
625  if( flag == infinitesimal ) then
626 
627  dstress(1:6) = matmul( d(1:6, 1:6), dstrain(1:6)-gauss%strain_bak(1:6) )
628  gauss%stress(1:6) = gauss%stress_bak(1:6) + dstress(1:6)
629  if( isviscoelastic(mtype) .AND. tincr /= 0.0d0 ) then
630  call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr, ttc, ttn, hdflag=hdflag_in )
631  gauss%stress = real(gauss%stress)
632  end if
633  gauss%strain_energy = gauss%strain_energy_bak+dot_product(gauss%stress_bak(1:6),dstrain(1:6)-gauss%strain_bak(1:6))
634  gauss%strain_energy = gauss%strain_energy+0.5d0*dot_product(dstress,dstrain(1:6)-gauss%strain_bak(1:6))
635 
636 
637  else if( flag == totallag ) then
638 
639  if( ishyperelastic(mtype) .OR. mtype == userelastic .OR. mtype == usermaterial ) then
640  call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, temp=0.d0, tempn=0.d0, hdflag=hdflag_in )
641  else if( ( isviscoelastic(mtype) .OR. mtype == norton ) .AND. tincr /= 0.0d0 ) then
642  gauss%pMaterial%mtype=mtype
643  call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr, ttc, ttn, hdflag=hdflag_in )
644  else
645  gauss%stress(1:6) = matmul( d(1:6, 1:6), dstrain(1:6) )
646  gauss%strain_energy = 0.5d0*dot_product(gauss%stress(1:6),dstrain(1:6))
647  end if
648 
649  else if( flag == updatelag ) then
650  ! CALL GEOMAT_C3( gauss%stress, mat )
651  ! D(:, :) = D(:, :)+mat(:, :)
652 
653  if( isviscoelastic(mtype) .AND. tincr /= 0.0d0 ) then
654  call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr, ttc, tt0, hdflag=hdflag_in )
655  else
656  dstress = real( matmul( d(1:6,1:6), dstrain(1:6) ) )
657  dumstress(1,1) = gauss%stress_bak(1)
658  dumstress(2,2) = gauss%stress_bak(2)
659  dumstress(3,3) = gauss%stress_bak(3)
660  dumstress(1,2) = gauss%stress_bak(4); dumstress(2,1)=dumstress(1,2)
661  dumstress(2,3) = gauss%stress_bak(5); dumstress(3,2)=dumstress(2,3)
662  dumstress(3,1) = gauss%stress_bak(6); dumstress(1,3)=dumstress(3,1)
663 
664  !stress integration
665  trd = dstrain(1)+dstrain(2)+dstrain(3)
666  dum(:,:) = dumstress + matmul( rot,dumstress ) - matmul( dumstress, rot ) - dumstress*trd
667  !call Hughes_Winget_rotation_3D( rot, dumstress, dum )
668 
669  gauss%stress(1) = dum(1,1) + dstress(1)
670  gauss%stress(2) = dum(2,2) + dstress(2)
671  gauss%stress(3) = dum(3,3) + dstress(3)
672  gauss%stress(4) = dum(1,2) + dstress(4)
673  gauss%stress(5) = dum(2,3) + dstress(5)
674  gauss%stress(6) = dum(3,1) + dstress(6)
675 
676  if( mtype == usermaterial ) then
677  call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, temp=0.d0, tempn=0.d0, hdflag=hdflag_in )
678  else if( mtype == norton ) then
679  !gauss%pMaterial%mtype = mtype
680  if( tincr /= 0.0d0 .AND. any( gauss%stress /= 0.0d0 ) ) then
681  !gauss%pMaterial%mtype = mtype
682  call stressupdate( gauss, d3, gauss%strain, gauss%stress, coordsys, time, tincr, ttc, ttn, hdflag=hdflag_in )
683  end if
684  end if
685  end if
686  gauss%strain_energy = gauss%strain_energy_bak+dot_product(gauss%stress(1:6),dstrain(1:6)) !s_t*de+de*C*de
687  gauss%strain_energy = gauss%strain_energy-0.5d0*dot_product(dstress(1:6),dstrain(1:6))
688 
689  end if
690 
691  if( iselastoplastic(mtype) ) then
692  call backwardeuler( gauss%pMaterial, gauss%stress, gauss%plstrain, &
693  gauss%istatus(1), gauss%fstatus, gauss%plpotential, ttc )
694  gauss%strain_energy = gauss%strain_energy+gauss%plpotential
695  end if
696 
697  !convert stress/strain measure for output
698  if( opsstype == kopss_solution ) then
699 
700  if( flag == infinitesimal ) then !linear
701  gauss%stress_out(1:6) = gauss%stress(1:6)
702  gauss%strain_out(1:6) = gauss%strain(1:6)
703  else !nonlinear
704  !convert stress
705  if( flag == totallag ) then
706  dumstress(1,1) = gauss%stress(1)
707  dumstress(2,2) = gauss%stress(2)
708  dumstress(3,3) = gauss%stress(3)
709  dumstress(1,2) = gauss%stress(4); dumstress(2,1)=dumstress(1,2)
710  dumstress(2,3) = gauss%stress(5); dumstress(3,2)=dumstress(2,3)
711  dumstress(3,1) = gauss%stress(6); dumstress(1,3)=dumstress(3,1)
712 
713  det = determinant33(f)
714  if( det == 0.d0 ) stop "Fail to convert stress: detF=0"
715  ! cauchy stress = (1/detF)*F*(2ndPK stress)*F^T
716  dumstress(1:3,1:3) = matmul(dumstress(1:3,1:3),transpose(f(1:3,1:3)))
717  dumstress(1:3,1:3) = (1.d0/det)*matmul(f(1:3,1:3),dumstress(1:3,1:3))
718 
719  gauss%stress_out(1) = dumstress(1,1)
720  gauss%stress_out(2) = dumstress(2,2)
721  gauss%stress_out(3) = dumstress(3,3)
722  gauss%stress_out(4) = dumstress(1,2)
723  gauss%stress_out(5) = dumstress(2,3)
724  gauss%stress_out(6) = dumstress(3,1)
725  else if( flag == updatelag ) then
726  gauss%stress_out(1:6) = gauss%stress(1:6)
727  endif
728 
729  !calc logarithmic strain
730  dum(1:3,1:3) = matmul(f(1:3,1:3),transpose(f(1:3,1:3)))
731  tensor(1) = dum(1,1)
732  tensor(2) = dum(2,2)
733  tensor(3) = dum(3,3)
734  tensor(4) = dum(1,2)
735  tensor(5) = dum(2,3)
736  tensor(6) = dum(3,1)
737  call get_principal(tensor, eigval, princ)
738 
739  do k=1,3
740  if( eigval(k) <= 0.d0 ) stop "Fail to calc log strain: stretch<0"
741  eigval(k) = 0.5d0*dlog(eigval(k)) !log(sqrt(lambda))
742  norm = dsqrt(dot_product(princ(1:3,k),princ(1:3,k)))
743  if( norm <= 0.d0 ) stop "Fail to calc log strain: stretch direction vector=0"
744  princ(1:3,k) = princ(1:3,k)/norm
745  end do
746  do i=1,3
747  do j=1,3
748  dum(i,j) = 0.d0
749  do k=1,3
750  dum(i,j) = dum(i,j) + eigval(k)*princ(i,k)*princ(j,k)
751  end do
752  end do
753  end do
754  gauss%strain_out(1) = dum(1,1)
755  gauss%strain_out(2) = dum(2,2)
756  gauss%strain_out(3) = dum(3,3)
757  gauss%strain_out(4) = 2.d0*dum(1,2)
758  gauss%strain_out(5) = 2.d0*dum(2,3)
759  gauss%strain_out(6) = 2.d0*dum(3,1)
760  endif
761 
762  else
763  gauss%stress_out(1:6) = gauss%stress(1:6)
764  gauss%strain_out(1:6) = gauss%strain(1:6)
765  end if
766 
767  end subroutine
768 
770  !---------------------------------------------------------------------*
771  subroutine update_c3 &
772  (etype,nn,ecoord, u, ddu, cdsys_id, coords, qf, &
773  gausses, iter, time, tincr, tt, t0, tn)
774  !---------------------------------------------------------------------*
775 
776  use m_fstr
777  use mmaterial
778  use mmechgauss
779  use m_matmatrix
780  use m_elastoplastic
781  use m_utilities
782 
783  integer(kind=kint), intent(in) :: etype
784  integer(kind=kint), intent(in) :: nn
785  real(kind=kreal), intent(in) :: ecoord(3, nn)
786  real(kind=kreal), intent(in) :: u(3, nn)
787  real(kind=kreal), intent(in) :: ddu(3, nn)
788  integer(kind=kint), intent(in) :: cdsys_id
789  real(kind=kreal), intent(inout) :: coords(3, 3)
790  real(kind=kreal), intent(out) :: qf(nn*3)
791  type(tgaussstatus), intent(inout) :: gausses(:)
792  integer, intent(in) :: iter
793  real(kind=kreal), intent(in) :: time
794  real(kind=kreal), intent(in) :: tincr
795  real(kind=kreal), intent(in) :: tt(nn)
796  real(kind=kreal), intent(in) :: t0(nn)
797  real(kind=kreal), intent(in) :: tn(nn)
798 
799  ! LOCAL VARIABLES
800  integer(kind=kint) :: flag
801  integer(kind=kint), parameter :: ndof = 3
802  real(kind=kreal) :: b(6,ndof*nn), b1(6,ndof*nn), spfunc(nn), ina(1)
803  real(kind=kreal) :: gderiv(nn,3), gderiv1(nn,3), gdispderiv(3,3), f(3,3), det, det1, wg, ttc, tt0, ttn
804  integer(kind=kint) :: j, lx, serr
805  real(kind=kreal) :: naturalcoord(3), rot(3,3), epsth(6)
806  real(kind=kreal) :: totaldisp(3,nn), elem(3,nn), elem1(3,nn), coordsys(3,3)
807  real(kind=kreal) :: dstrain(6)
808  real(kind=kreal) :: alpo(3)
809  logical :: ierr, matlaniso
810 
811  qf(:) = 0.0d0
812  ! we suppose the same material type in the element
813  flag = gausses(1)%pMaterial%nlgeom_flag
814  elem(:,:) = ecoord(:,:)
815  totaldisp(:,:) = u(:,:)+ ddu(:,:)
816  if( flag == updatelag ) then
817  elem(:,:) = (0.5d0*ddu(:,:)+u(:,:) ) +ecoord(:,:)
818  elem1(:,:) = (ddu(:,:)+u(:,:) ) +ecoord(:,:)
819  ! elem = elem1
820  totaldisp(:,:) = ddu(:,:)
821  end if
822 
823  matlaniso = .false.
824  ina = tt(1)
825  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
826  if( .not. ierr ) matlaniso = .true.
827 
828  do lx = 1, numofquadpoints(etype)
829 
830  call getquadpoint( etype, lx, naturalcoord(:) )
831  call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
832 
833  if( cdsys_id > 0 ) then
834  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:,:), serr )
835  if( serr == -1 ) stop "Fail to setup local coordinate"
836  if( serr == -2 ) then
837  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
838  end if
839  end if
840 
841  ! ========================================================
842  ! UPDATE STRAIN and STRESS
843  ! ========================================================
844 
845  ! Thermal Strain
846  epsth = 0.0d0
847  call getshapefunc(etype, naturalcoord, spfunc)
848  ttc = dot_product(tt, spfunc)
849  tt0 = dot_product(t0, spfunc)
850  ttn = dot_product(tn, spfunc)
851  call cal_thermal_expansion_c3( tt0, ttc, gausses(lx)%pMaterial, coordsys, matlaniso, epsth )
852 
853  ! Update strain
854  ! Small strain
855  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
856  dstrain(1) = gdispderiv(1, 1)
857  dstrain(2) = gdispderiv(2, 2)
858  dstrain(3) = gdispderiv(3, 3)
859  dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
860  dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
861  dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
862  dstrain(:) = dstrain(:)-epsth(:) ! alright?
863 
864  f(1:3,1:3) = 0.d0; f(1,1)=1.d0; f(2,2)=1.d0; f(3,3)=1.d0; !deformation gradient
865  if( flag == infinitesimal ) then
866  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(1:6)
867 
868  else if( flag == totallag ) then
869  ! Green-Lagrange strain
870  dstrain(1) = dstrain(1)+0.5d0*dot_product( gdispderiv(:, 1), gdispderiv(:, 1) )
871  dstrain(2) = dstrain(2)+0.5d0*dot_product( gdispderiv(:, 2), gdispderiv(:, 2) )
872  dstrain(3) = dstrain(3)+0.5d0*dot_product( gdispderiv(:, 3), gdispderiv(:, 3) )
873  dstrain(4) = dstrain(4)+( gdispderiv(1, 1)*gdispderiv(1, 2) &
874  +gdispderiv(2, 1)*gdispderiv(2, 2)+gdispderiv(3, 1)*gdispderiv(3, 2) )
875  dstrain(5) = dstrain(5)+( gdispderiv(1, 2)*gdispderiv(1, 3) &
876  +gdispderiv(2, 2)*gdispderiv(2, 3)+gdispderiv(3, 2)*gdispderiv(3, 3) )
877  dstrain(6) = dstrain(6)+( gdispderiv(1, 1)*gdispderiv(1, 3) &
878  +gdispderiv(2, 1)*gdispderiv(2, 3)+gdispderiv(3, 1)*gdispderiv(3, 3) )
879 
880  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
881  f(1:3,1:3) = f(1:3,1:3) + gdispderiv(1:3,1:3)
882 
883  else if( flag == updatelag ) then
884  rot = 0.0d0
885  rot(1, 2)= 0.5d0*(gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
886  rot(2, 3)= 0.5d0*(gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
887  rot(1, 3)= 0.5d0*(gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
888 
889  gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
890 
891  call getglobalderiv(etype, nn, naturalcoord, ecoord, det1, gderiv1)
892  f(1:3,1:3) = f(1:3,1:3) + matmul( u(1:ndof, 1:nn)+ddu(1:ndof, 1:nn), gderiv1(1:nn, 1:ndof) )
893 
894  end if
895 
896  ! Update stress
897  call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn )
898 
899  ! ========================================================
900  ! calculate the internal force ( equivalent nodal force )
901  ! ========================================================
902  ! Small strain
903  b(1:6, 1:nn*ndof) = 0.0d0
904  do j=1,nn
905  b(1,3*j-2) = gderiv(j, 1)
906  b(2,3*j-1) = gderiv(j, 2)
907  b(3,3*j ) = gderiv(j, 3)
908  b(4,3*j-2) = gderiv(j, 2)
909  b(4,3*j-1) = gderiv(j, 1)
910  b(5,3*j-1) = gderiv(j, 3)
911  b(5,3*j ) = gderiv(j, 2)
912  b(6,3*j-2) = gderiv(j, 3)
913  b(6,3*j ) = gderiv(j, 1)
914  end do
915 
916  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
917  if( flag == infinitesimal ) then
918 
919  else if( flag == totallag ) then
920 
921  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
922  b1(1:6, 1:nn*ndof)=0.0d0
923  do j = 1,nn
924  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
925  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
926  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
927  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
928  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
929  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
930  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
931  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
932  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
933  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
934  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
935  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
936  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
937  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
938  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
939  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
940  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
941  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
942  end do
943  ! BL = BL0 + BL1
944  do j=1,nn*ndof
945  b(:,j) = b(:,j)+b1(:,j)
946  end do
947 
948  else if( flag == updatelag ) then
949 
950  call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv)
951  b(1:6, 1:nn*ndof) = 0.0d0
952  do j = 1, nn
953  b(1, 3*j-2) = gderiv(j, 1)
954  b(2, 3*j-1) = gderiv(j, 2)
955  b(3, 3*j ) = gderiv(j, 3)
956  b(4, 3*j-2) = gderiv(j, 2)
957  b(4, 3*j-1) = gderiv(j, 1)
958  b(5, 3*j-1) = gderiv(j, 3)
959  b(5, 3*j ) = gderiv(j, 2)
960  b(6, 3*j-2) = gderiv(j, 3)
961  b(6, 3*j ) = gderiv(j, 1)
962  end do
963 
964  end if
965 
966  ! calculate the Internal Force
967  wg=getweight( etype, lx )*det
968  qf(1:nn*ndof) &
969  = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
970 
971  ! integrate strain energy
972  gausses(lx)%strain_energy = gausses(lx)%strain_energy*wg
973  end do
974 
975  end subroutine update_c3
976  !
977  !----------------------------------------------------------------------*
978  subroutine nodalstress_c3(etype, nn, gausses, ndstrain, ndstress)
979  !----------------------------------------------------------------------*
980  !
981  ! Calculate Strain and Stress increment of solid elements
982  !
983  use mmechgauss
984 
985  !---------------------------------------------------------------------
986 
987  integer(kind=kint), intent(in) :: etype, nn
988  type(tgaussstatus), intent(in) :: gausses(:)
989  real(kind=kreal), intent(out) :: ndstrain(nn,6)
990  real(kind=kreal), intent(out) :: ndstress(nn,6)
991 
992  !---------------------------------------------------------------------
993 
994  integer :: i, ic
995  real(kind=kreal) :: temp(12)
996 
997  !---------------------------------------------------------------------
998 
999  temp(:) = 0.0d0
1000 
1001  ic = numofquadpoints(etype)
1002 
1003  do i = 1, ic
1004  temp(1:6) = temp(1:6) +gausses(i)%strain_out(1:6)
1005  temp(7:12) = temp(7:12)+gausses(i)%stress_out(1:6)
1006  end do
1007 
1008  temp(1:12) = temp(1:12)/ic
1009 
1010  do i=1,nn
1011  ndstrain(i, 1:6) = temp(1:6)
1012  ndstress(i, 1:6) = temp(7:12)
1013  enddo
1014 
1015  end subroutine nodalstress_c3
1016 
1017 
1018  !----------------------------------------------------------------------*
1019  subroutine elementstress_c3(etype, gausses, strain, stress)
1020  !----------------------------------------------------------------------*
1021  !
1022  ! Calculate Strain and Stress increment of solid elements
1023  !
1024  use mmechgauss
1025 
1026  !---------------------------------------------------------------------
1027 
1028  integer(kind=kint), intent(in) :: etype
1029  type(tgaussstatus), intent(in) :: gausses(:)
1030  real(kind=kreal), intent(out) :: strain(6)
1031  real(kind=kreal), intent(out) :: stress(6)
1032 
1033  !---------------------------------------------------------------------
1034 
1035  integer :: i, ic
1036 
1037  !---------------------------------------------------------------------
1038 
1039  strain(:) = 0.0d0; stress(:) = 0.0d0
1040 
1041  ic = numofquadpoints(etype)
1042 
1043  do i = 1, ic
1044  strain(:) = strain(:)+gausses(i)%strain_out(1:6)
1045  stress(:) = stress(:)+gausses(i)%stress_out(1:6)
1046  enddo
1047 
1048  strain(:) = strain(:)/ic
1049  stress(:) = stress(:)/ic
1050 
1051  end subroutine elementstress_c3
1052 
1053 
1055  !----------------------------------------------------------------------*
1056  real(kind=kreal) function volume_c3(etype, nn, XX, YY, ZZ)
1057  !----------------------------------------------------------------------*
1058 
1059  integer(kind=kint), intent(in) :: etype, nn
1060  real(kind=kreal), intent(in) :: xx(:), yy(:), zz(:)
1061 
1062  !---------------------------------------------------------------------
1063 
1064  real(kind=kreal) :: xj(3, 3), det
1065  integer(kind=kint) :: lx
1066  real(kind=kreal) :: localcoord(3), deriv(nn, 3)
1067 
1068  !---------------------------------------------------------------------
1069 
1070  volume_c3 = 0.0d0
1071 
1072  ! LOOP FOR INTEGRATION POINTS
1073  do lx = 1, numofquadpoints(etype)
1074 
1075  call getquadpoint(etype, lx, localcoord)
1076  call getshapederiv(etype, localcoord, deriv)
1077 
1078  ! JACOBI MATRIX
1079  xj(1, 1:3)= matmul( xx(1:nn), deriv(1:nn,1:3) )
1080  xj(2, 1:3)= matmul( yy(1:nn), deriv(1:nn,1:3) )
1081  xj(3, 1:3)= matmul( zz(1:nn), deriv(1:nn,1:3) )
1082 
1083  ! DETERMINANT OF JACOBIAN
1084  det = xj(1, 1)*xj(2, 2)*xj(3, 3) &
1085  +xj(2, 1)*xj(3, 2)*xj(1, 3) &
1086  +xj(3, 1)*xj(1, 2)*xj(2, 3) &
1087  -xj(3, 1)*xj(2, 2)*xj(1, 3) &
1088  -xj(2, 1)*xj(1, 2)*xj(3, 3) &
1089  -xj(1, 1)*xj(3, 2)*xj(2, 3)
1090 
1091  volume_c3 = volume_c3+getweight(etype, lx)*det
1092 
1093  end do
1094 
1095  end function volume_c3
1096 
1097 
1098 end module m_static_lib_3d
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 coordinate system.
Definition: element.f90:647
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
Definition: element.f90:131
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:489
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:741
subroutine getsubface(intype, innumber, outtype, nodes)
Find the definition of surface of the element.
Definition: element.f90:193
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate derivatives of shape function in natural coordinate system.
Definition: element.f90:578
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:535
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:450
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
Definition: element.f90:870
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 defines common data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint), parameter kopss_solution
Definition: m_fstr.f90:130
integer(kind=kint) opsstype
Definition: m_fstr.f90:132
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.f90:136
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
subroutine matlmatrix(gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp, hdflag)
Calculate constituive matrix.
subroutine stressupdate(gauss, sectType, strain, stress, cdsys, time, dtime, temp, tempn, hdflag)
Update strain and stress for elastic and hyperelastic materials.
This module provides common functions of Solid elements.
subroutine elementstress_c3(etype, gausses, strain, stress)
subroutine tload_c3(etype, nn, XX, YY, ZZ, TT, T0, gausses, VECT, cdsys_ID, coords)
This subroutien calculate thermal loading.
subroutine dl_c3(etype, nn, XX, YY, ZZ, RHO, LTYPE, PARAMS, VECT, NSIZE)
Distributed external load.
subroutine hughes_winget_rotation_3d(rot, stress_in, stress_out)
subroutine geomat_c3(stress, mat)
subroutine update_c3(etype, nn, ecoord, u, ddu, cdsys_ID, coords, qf, gausses, iter, time, tincr, TT, T0, TN)
Update strain and stress inside element.
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, hdflag)
subroutine stf_c3(etype, nn, ecoord, gausses, stiff, cdsys_ID, coords, time, tincr, u, temperature)
This subroutine calculate stiff matrix of general solid elements.
real(kind=kreal) function volume_c3(etype, nn, XX, YY, ZZ)
Volume of element.
subroutine nodalstress_c3(etype, nn, gausses, ndstrain, ndstress)
This module provides aux functions.
Definition: utilities.f90:6
real(kind=kreal) function determinant33(XJ)
Compute determinant for 3*3 matrix.
Definition: utilities.f90:270
subroutine get_principal(tensor, eigval, princmatrix)
Definition: utilities.f90:411
subroutine transformation(jacob, tm)
Definition: utilities.f90:376
subroutine calinverse(NN, A)
calculate inverse of matrix a
Definition: utilities.f90:295
This module summarizes all information of material properties.
Definition: material.f90:6
integer(kind=kint), parameter totallag
Definition: material.f90:14
character(len=dict_key_length) mc_orthoexp
Definition: material.f90:144
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