FrontISTR  5.9.0
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 !-------------------------------------------------------------------------------
6 module m_static_lib_3d
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, log_valid
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  if( any(eigval(1:3) <= 0.d0) ) then
740  gauss%strain_out(1:6) = gauss%strain(1:6)
741  else
742  log_valid = 1
743  do k=1,3
744  norm = dsqrt(dot_product(princ(1:3,k),princ(1:3,k)))
745  if( norm <= 1.0d-15 ) then
746  log_valid = 0
747  exit
748  end if
749  eigval(k) = 0.5d0*dlog(eigval(k))
750  princ(1:3,k) = princ(1:3,k)/norm
751  end do
752  if( log_valid == 0 ) then
753  gauss%strain_out(1:6) = gauss%strain(1:6)
754  else
755  do i=1,3
756  do j=1,3
757  dum(i,j) = 0.d0
758  do k=1,3
759  dum(i,j) = dum(i,j) + eigval(k)*princ(i,k)*princ(j,k)
760  end do
761  end do
762  end do
763  gauss%strain_out(1) = dum(1,1)
764  gauss%strain_out(2) = dum(2,2)
765  gauss%strain_out(3) = dum(3,3)
766  gauss%strain_out(4) = 2.d0*dum(1,2)
767  gauss%strain_out(5) = 2.d0*dum(2,3)
768  gauss%strain_out(6) = 2.d0*dum(3,1)
769  end if
770  end if
771  endif
772 
773  else
774  gauss%stress_out(1:6) = gauss%stress(1:6)
775  gauss%strain_out(1:6) = gauss%strain(1:6)
776  end if
777 
778  end subroutine
779 
781  !---------------------------------------------------------------------*
782  subroutine update_c3 &
783  (etype,nn,ecoord, u, ddu, cdsys_id, coords, qf, &
784  gausses, iter, time, tincr, tt, t0, tn)
785  !---------------------------------------------------------------------*
786 
787  use m_fstr
788  use mmaterial
789  use mmechgauss
790  use m_matmatrix
791  use m_elastoplastic
792  use m_utilities
793 
794  integer(kind=kint), intent(in) :: etype
795  integer(kind=kint), intent(in) :: nn
796  real(kind=kreal), intent(in) :: ecoord(3, nn)
797  real(kind=kreal), intent(in) :: u(3, nn)
798  real(kind=kreal), intent(in) :: ddu(3, nn)
799  integer(kind=kint), intent(in) :: cdsys_id
800  real(kind=kreal), intent(inout) :: coords(3, 3)
801  real(kind=kreal), intent(out) :: qf(nn*3)
802  type(tgaussstatus), intent(inout) :: gausses(:)
803  integer, intent(in) :: iter
804  real(kind=kreal), intent(in) :: time
805  real(kind=kreal), intent(in) :: tincr
806  real(kind=kreal), intent(in) :: tt(nn)
807  real(kind=kreal), intent(in) :: t0(nn)
808  real(kind=kreal), intent(in) :: tn(nn)
809 
810  ! LOCAL VARIABLES
811  integer(kind=kint) :: flag
812  integer(kind=kint), parameter :: ndof = 3
813  real(kind=kreal) :: b(6,ndof*nn), b1(6,ndof*nn), spfunc(nn), ina(1)
814  real(kind=kreal) :: gderiv(nn,3), gderiv1(nn,3), gdispderiv(3,3), f(3,3), det, det1, wg, ttc, tt0, ttn
815  integer(kind=kint) :: j, lx, serr
816  real(kind=kreal) :: naturalcoord(3), rot(3,3), epsth(6)
817  real(kind=kreal) :: totaldisp(3,nn), elem(3,nn), elem1(3,nn), coordsys(3,3)
818  real(kind=kreal) :: dstrain(6)
819  real(kind=kreal) :: alpo(3)
820  logical :: ierr, matlaniso
821 
822  qf(:) = 0.0d0
823  ! we suppose the same material type in the element
824  flag = gausses(1)%pMaterial%nlgeom_flag
825  elem(:,:) = ecoord(:,:)
826  totaldisp(:,:) = u(:,:)+ ddu(:,:)
827  if( flag == updatelag ) then
828  elem(:,:) = (0.5d0*ddu(:,:)+u(:,:) ) +ecoord(:,:)
829  elem1(:,:) = (ddu(:,:)+u(:,:) ) +ecoord(:,:)
830  ! elem = elem1
831  totaldisp(:,:) = ddu(:,:)
832  end if
833 
834  matlaniso = .false.
835  ina = tt(1)
836  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
837  if( .not. ierr ) matlaniso = .true.
838 
839  do lx = 1, numofquadpoints(etype)
840 
841  call getquadpoint( etype, lx, naturalcoord(:) )
842  call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
843 
844  if( cdsys_id > 0 ) then
845  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:,:), serr )
846  if( serr == -1 ) stop "Fail to setup local coordinate"
847  if( serr == -2 ) then
848  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
849  end if
850  end if
851 
852  ! ========================================================
853  ! UPDATE STRAIN and STRESS
854  ! ========================================================
855 
856  ! Thermal Strain
857  epsth = 0.0d0
858  call getshapefunc(etype, naturalcoord, spfunc)
859  ttc = dot_product(tt, spfunc)
860  tt0 = dot_product(t0, spfunc)
861  ttn = dot_product(tn, spfunc)
862  call cal_thermal_expansion_c3( tt0, ttc, gausses(lx)%pMaterial, coordsys, matlaniso, epsth )
863 
864  ! Update strain
865  ! Small strain
866  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
867  dstrain(1) = gdispderiv(1, 1)
868  dstrain(2) = gdispderiv(2, 2)
869  dstrain(3) = gdispderiv(3, 3)
870  dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
871  dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
872  dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
873  dstrain(:) = dstrain(:)-epsth(:) ! alright?
874 
875  f(1:3,1:3) = 0.d0; f(1,1)=1.d0; f(2,2)=1.d0; f(3,3)=1.d0; !deformation gradient
876  if( flag == infinitesimal ) then
877  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(1:6)
878 
879  else if( flag == totallag ) then
880  ! Green-Lagrange strain
881  dstrain(1) = dstrain(1)+0.5d0*dot_product( gdispderiv(:, 1), gdispderiv(:, 1) )
882  dstrain(2) = dstrain(2)+0.5d0*dot_product( gdispderiv(:, 2), gdispderiv(:, 2) )
883  dstrain(3) = dstrain(3)+0.5d0*dot_product( gdispderiv(:, 3), gdispderiv(:, 3) )
884  dstrain(4) = dstrain(4)+( gdispderiv(1, 1)*gdispderiv(1, 2) &
885  +gdispderiv(2, 1)*gdispderiv(2, 2)+gdispderiv(3, 1)*gdispderiv(3, 2) )
886  dstrain(5) = dstrain(5)+( gdispderiv(1, 2)*gdispderiv(1, 3) &
887  +gdispderiv(2, 2)*gdispderiv(2, 3)+gdispderiv(3, 2)*gdispderiv(3, 3) )
888  dstrain(6) = dstrain(6)+( gdispderiv(1, 1)*gdispderiv(1, 3) &
889  +gdispderiv(2, 1)*gdispderiv(2, 3)+gdispderiv(3, 1)*gdispderiv(3, 3) )
890 
891  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
892  f(1:3,1:3) = f(1:3,1:3) + gdispderiv(1:3,1:3)
893 
894  else if( flag == updatelag ) then
895  rot = 0.0d0
896  rot(1, 2)= 0.5d0*(gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
897  rot(2, 3)= 0.5d0*(gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
898  rot(1, 3)= 0.5d0*(gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
899 
900  gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
901 
902  call getglobalderiv(etype, nn, naturalcoord, ecoord, det1, gderiv1)
903  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) )
904 
905  end if
906 
907  ! Update stress
908  call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn )
909 
910  ! ========================================================
911  ! calculate the internal force ( equivalent nodal force )
912  ! ========================================================
913  ! Small strain
914  b(1:6, 1:nn*ndof) = 0.0d0
915  do j=1,nn
916  b(1,3*j-2) = gderiv(j, 1)
917  b(2,3*j-1) = gderiv(j, 2)
918  b(3,3*j ) = gderiv(j, 3)
919  b(4,3*j-2) = gderiv(j, 2)
920  b(4,3*j-1) = gderiv(j, 1)
921  b(5,3*j-1) = gderiv(j, 3)
922  b(5,3*j ) = gderiv(j, 2)
923  b(6,3*j-2) = gderiv(j, 3)
924  b(6,3*j ) = gderiv(j, 1)
925  end do
926 
927  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
928  if( flag == infinitesimal ) then
929 
930  else if( flag == totallag ) then
931 
932  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
933  b1(1:6, 1:nn*ndof)=0.0d0
934  do j = 1,nn
935  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
936  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
937  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
938  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
939  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
940  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
941  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
942  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
943  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
944  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
945  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
946  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
947  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
948  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
949  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
950  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
951  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
952  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
953  end do
954  ! BL = BL0 + BL1
955  do j=1,nn*ndof
956  b(:,j) = b(:,j)+b1(:,j)
957  end do
958 
959  else if( flag == updatelag ) then
960 
961  call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv)
962  b(1:6, 1:nn*ndof) = 0.0d0
963  do j = 1, nn
964  b(1, 3*j-2) = gderiv(j, 1)
965  b(2, 3*j-1) = gderiv(j, 2)
966  b(3, 3*j ) = gderiv(j, 3)
967  b(4, 3*j-2) = gderiv(j, 2)
968  b(4, 3*j-1) = gderiv(j, 1)
969  b(5, 3*j-1) = gderiv(j, 3)
970  b(5, 3*j ) = gderiv(j, 2)
971  b(6, 3*j-2) = gderiv(j, 3)
972  b(6, 3*j ) = gderiv(j, 1)
973  end do
974 
975  end if
976 
977  ! calculate the Internal Force
978  wg=getweight( etype, lx )*det
979  qf(1:nn*ndof) &
980  = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
981 
982  ! integrate strain energy
983  gausses(lx)%strain_energy = gausses(lx)%strain_energy*wg
984  end do
985 
986  end subroutine update_c3
987  !
988  !----------------------------------------------------------------------*
989  subroutine nodalstress_c3(etype, nn, gausses, ndstrain, ndstress)
990  !----------------------------------------------------------------------*
991  !
992  ! Calculate Strain and Stress increment of solid elements
993  !
994  use mmechgauss
995 
996  !---------------------------------------------------------------------
997 
998  integer(kind=kint), intent(in) :: etype, nn
999  type(tgaussstatus), intent(in) :: gausses(:)
1000  real(kind=kreal), intent(out) :: ndstrain(nn,6)
1001  real(kind=kreal), intent(out) :: ndstress(nn,6)
1002 
1003  !---------------------------------------------------------------------
1004 
1005  integer :: i, ic
1006  real(kind=kreal) :: temp(12)
1007 
1008  !---------------------------------------------------------------------
1009 
1010  temp(:) = 0.0d0
1011 
1012  ic = numofquadpoints(etype)
1013 
1014  do i = 1, ic
1015  temp(1:6) = temp(1:6) +gausses(i)%strain_out(1:6)
1016  temp(7:12) = temp(7:12)+gausses(i)%stress_out(1:6)
1017  end do
1018 
1019  temp(1:12) = temp(1:12)/ic
1020 
1021  do i=1,nn
1022  ndstrain(i, 1:6) = temp(1:6)
1023  ndstress(i, 1:6) = temp(7:12)
1024  enddo
1025 
1026  end subroutine nodalstress_c3
1027 
1028 
1029  !----------------------------------------------------------------------*
1030  subroutine elementstress_c3(etype, gausses, strain, stress)
1031  !----------------------------------------------------------------------*
1032  !
1033  ! Calculate Strain and Stress increment of solid elements
1034  !
1035  use mmechgauss
1036 
1037  !---------------------------------------------------------------------
1038 
1039  integer(kind=kint), intent(in) :: etype
1040  type(tgaussstatus), intent(in) :: gausses(:)
1041  real(kind=kreal), intent(out) :: strain(6)
1042  real(kind=kreal), intent(out) :: stress(6)
1043 
1044  !---------------------------------------------------------------------
1045 
1046  integer :: i, ic
1047 
1048  !---------------------------------------------------------------------
1049 
1050  strain(:) = 0.0d0; stress(:) = 0.0d0
1051 
1052  ic = numofquadpoints(etype)
1053 
1054  do i = 1, ic
1055  strain(:) = strain(:)+gausses(i)%strain_out(1:6)
1056  stress(:) = stress(:)+gausses(i)%stress_out(1:6)
1057  enddo
1058 
1059  strain(:) = strain(:)/ic
1060  stress(:) = stress(:)/ic
1061 
1062  end subroutine elementstress_c3
1063 
1064 
1066  !----------------------------------------------------------------------*
1067  real(kind=kreal) function volume_c3(etype, nn, XX, YY, ZZ)
1068  !----------------------------------------------------------------------*
1069 
1070  integer(kind=kint), intent(in) :: etype, nn
1071  real(kind=kreal), intent(in) :: xx(:), yy(:), zz(:)
1072 
1073  !---------------------------------------------------------------------
1074 
1075  real(kind=kreal) :: xj(3, 3), det
1076  integer(kind=kint) :: lx
1077  real(kind=kreal) :: localcoord(3), deriv(nn, 3)
1078 
1079  !---------------------------------------------------------------------
1080 
1081  volume_c3 = 0.0d0
1082 
1083  ! LOOP FOR INTEGRATION POINTS
1084  do lx = 1, numofquadpoints(etype)
1085 
1086  call getquadpoint(etype, lx, localcoord)
1087  call getshapederiv(etype, localcoord, deriv)
1088 
1089  ! JACOBI MATRIX
1090  xj(1, 1:3)= matmul( xx(1:nn), deriv(1:nn,1:3) )
1091  xj(2, 1:3)= matmul( yy(1:nn), deriv(1:nn,1:3) )
1092  xj(3, 1:3)= matmul( zz(1:nn), deriv(1:nn,1:3) )
1093 
1094  ! DETERMINANT OF JACOBIAN
1095  det = xj(1, 1)*xj(2, 2)*xj(3, 3) &
1096  +xj(2, 1)*xj(3, 2)*xj(1, 3) &
1097  +xj(3, 1)*xj(1, 2)*xj(2, 3) &
1098  -xj(3, 1)*xj(2, 2)*xj(1, 3) &
1099  -xj(2, 1)*xj(1, 2)*xj(3, 3) &
1100  -xj(1, 1)*xj(3, 2)*xj(2, 3)
1101 
1102  volume_c3 = volume_c3+getweight(etype, lx)*det
1103 
1104  end do
1105 
1106  end function volume_c3
1107 
1108 
1109 end module m_static_lib_3d
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::getnumberofnodes
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
Definition: element.f90:131
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::getsubface
subroutine getsubface(intype, innumber, outtype, nodes)
Find the definition of surface of the element.
Definition: element.f90:193
elementinfo::getshapederiv
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate derivatives of shape function in natural coordinate system.
Definition: element.f90:578
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
elementinfo::surfacenormal
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
Definition: element.f90:870
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_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.F90:15
m_fstr::kopss_solution
integer(kind=kint), parameter kopss_solution
Definition: m_fstr.F90:132
m_fstr::opsstype
integer(kind=kint) opsstype
Definition: m_fstr.F90:134
m_fstr::ref_temp
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.F90:138
m_matmatrix
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
m_matmatrix::matlmatrix
subroutine matlmatrix(gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp, hdflag)
Calculate constituive matrix.
Definition: calMatMatrix.f90:30
m_matmatrix::stressupdate
subroutine stressupdate(gauss, sectType, strain, stress, cdsys, time, dtime, temp, tempn, hdflag)
Update strain and stress for elastic and hyperelastic materials.
Definition: calMatMatrix.f90:103
m_static_lib_3d
This module provides common functions of Solid elements.
Definition: static_LIB_3d.f90:6
m_static_lib_3d::elementstress_c3
subroutine elementstress_c3(etype, gausses, strain, stress)
Definition: static_LIB_3d.f90:1031
m_static_lib_3d::tload_c3
subroutine tload_c3(etype, nn, XX, YY, ZZ, TT, T0, gausses, VECT, cdsys_ID, coords)
This subroutien calculate thermal loading.
Definition: static_LIB_3d.f90:389
m_static_lib_3d::dl_c3
subroutine dl_c3(etype, nn, XX, YY, ZZ, RHO, LTYPE, PARAMS, VECT, NSIZE)
Distributed external load.
Definition: static_LIB_3d.f90:214
m_static_lib_3d::hughes_winget_rotation_3d
subroutine hughes_winget_rotation_3d(rot, stress_in, stress_out)
Definition: static_LIB_3d.f90:570
m_static_lib_3d::geomat_c3
subroutine geomat_c3(stress, mat)
Definition: static_LIB_3d.f90:19
m_static_lib_3d::update_c3
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.
Definition: static_LIB_3d.f90:785
m_static_lib_3d::cal_thermal_expansion_c3
subroutine cal_thermal_expansion_c3(tt0, ttc, material, coordsys, matlaniso, EPSTH)
Definition: static_LIB_3d.f90:519
m_static_lib_3d::update_stress3d
subroutine update_stress3d(flag, gauss, rot, dstrain, F, coordsys, time, tincr, ttc, tt0, ttn, hdflag)
Definition: static_LIB_3d.f90:587
m_static_lib_3d::stf_c3
subroutine stf_c3(etype, nn, ecoord, gausses, stiff, cdsys_ID, coords, time, tincr, u, temperature)
This subroutine calculate stiff matrix of general solid elements.
Definition: static_LIB_3d.f90:53
m_static_lib_3d::volume_c3
real(kind=kreal) function volume_c3(etype, nn, XX, YY, ZZ)
Volume of element.
Definition: static_LIB_3d.f90:1068
m_static_lib_3d::nodalstress_c3
subroutine nodalstress_c3(etype, nn, gausses, ndstrain, ndstress)
Definition: static_LIB_3d.f90:990
m_utilities
This module provides aux functions.
Definition: utilities.f90:6
m_utilities::determinant33
real(kind=kreal) function determinant33(XJ)
Compute determinant for 3*3 matrix.
Definition: utilities.f90:270
m_utilities::get_principal
subroutine get_principal(tensor, eigval, princmatrix)
Definition: utilities.f90:450
m_utilities::transformation
subroutine transformation(jacob, tm)
Definition: utilities.f90:415
m_utilities::calinverse
subroutine calinverse(NN, A)
calculate inverse of matrix a
Definition: utilities.f90:295
mmaterial
This module summarizes all information of material properties.
Definition: material.f90:6
mmaterial::totallag
integer(kind=kint), parameter totallag
Definition: material.f90:14
mmaterial::mc_orthoexp
character(len=dict_key_length) mc_orthoexp
Definition: material.f90:144
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:16