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