FrontISTR  5.2.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 !-------------------------------------------------------------------------------
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), optional :: 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  if( present(temperature) ) then
110  call getshapefunc(etype, naturalcoord, spfunc)
111  temp = dot_product(temperature, spfunc)
112  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys, temp )
113  else
114  call matlmatrix( gausses(lx), d3, d, time, tincr, coordsys )
115  end if
116 
117  if( flag == updatelag ) then
118  call geomat_c3( gausses(lx)%stress, mat )
119  d(:, :) = d(:, :)-mat
120  endif
121 
122  wg = getweight(etype, lx)*det
123  b(1:6, 1:nn*ndof) = 0.0d0
124  do j = 1, nn
125  b(1, 3*j-2)=gderiv(j, 1)
126  b(2, 3*j-1)=gderiv(j, 2)
127  b(3, 3*j )=gderiv(j, 3)
128  b(4, 3*j-2)=gderiv(j, 2)
129  b(4, 3*j-1)=gderiv(j, 1)
130  b(5, 3*j-1)=gderiv(j, 3)
131  b(5, 3*j )=gderiv(j, 2)
132  b(6, 3*j-2)=gderiv(j, 3)
133  b(6, 3*j )=gderiv(j, 1)
134  enddo
135 
136  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
137  if( flag == totallag ) then
138  ! ---dudx(i, j) ==> gdispderiv(i, j)
139  gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
140  b1(1:6, 1:nn*ndof)=0.0d0
141  do j=1, nn
142  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
143  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
144  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
145  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
146  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
147  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
148  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
149  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
150  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
151  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
152  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
153  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
154  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
155  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
156  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
157  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
158  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
159  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
160  end do
161  ! ---BL = BL0 + BL1
162  do j=1, nn*ndof
163  b(:, j) = b(:, j)+b1(:, j)
164  end do
165  end if
166 
167  db(1:6, 1:nn*ndof) = matmul( d, b(1:6, 1:nn*ndof) )
168  forall( i=1:nn*ndof, j=1:nn*ndof )
169  stiff(i, j) = stiff(i, j)+dot_product( b(:, i), db(:, j) )*wg
170  end forall
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  forall( i=1:nn*ndof, j=1:nn*ndof )
201  stiff(i, j) = stiff(i, j)+dot_product( bn(:, i), sbn(:, j) )*wg
202  end forall
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  ! SELCTION 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*val*phx
338  coefy=rho*val*val*phy
339  coefz=rho*val*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 exapansion 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( matlaniso ) 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  ina(1) = ttc
534  if( matlaniso ) then
535  call fetch_tabledata( mc_orthoexp, material%dict, alpo(:), ierr, ina )
536  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
537  else
538  call fetch_tabledata( mc_themoexp, material%dict, outa(:), ierr, ina )
539  if( ierr ) outa(1) = material%variables(m_exapnsion)
540  alp = outa(1)
541  end if
542  ina(1) = tt0
543  if( matlaniso ) then
544  call fetch_tabledata( mc_orthoexp, material%dict, alpo0(:), ierr, ina )
545  if( ierr ) stop "Fails in fetching orthotropic expansion coefficient!"
546  else
547  call fetch_tabledata( mc_themoexp, material%dict, outa(:), ierr, ina )
548  if( ierr ) outa(1) = material%variables(m_exapnsion)
549  alp0 = outa(1)
550  end if
551  if( matlaniso ) then
552  do j=1,3
553  epsth(j) = alpo(j)*(ttc-ref_temp)-alpo0(j)*(tt0-ref_temp)
554  end do
555  call transformation( coordsys(:, :), tm)
556  epsth(:) = matmul( epsth(:), tm ) ! to global coord
557  epsth(4:6) = epsth(4:6)*2.0d0
558  else
559  epsth(1:3)=alp*(ttc-ref_temp)-alp0*(tt0-ref_temp)
560  end if
561  end subroutine
562 
563  !Hughes, T. J. R., and J. Winget, ?gFinite Rotation Effects in Numerical Integration of
564  ! Rate Constitutive Equations Arising in Large Deformation Analysis,?h
565  ! International Journal for Numerical Methods in Engineering, vol. 15, pp. 1862-1867, 1980.
566  subroutine hughes_winget_rotation_3d( rot, stress_in, stress_out )
567  use m_utilities
568  real(kind=kreal), intent(in) :: rot(3,3)
569  real(kind=kreal), intent(in) :: stress_in(3,3)
570  real(kind=kreal), intent(out) :: stress_out(3,3)
571 
572  real(kind=kreal) :: dr(3,3)
573 
574  !calc dR=(I-0.5*dW)^-1*(I+0.5*dW)
575  dr = i3-0.5d0*rot
576  call calinverse(3, dr)
577  dr = matmul(dr,i3+0.5d0*rot)
578 
579  stress_out = matmul(dr,stress_in)
580  stress_out = matmul(stress_out,transpose(dr))
581  end subroutine
582 
583  subroutine update_stress3d( flag, gauss, rot, dstrain, F, coordsys, time, tincr, ttc, tt0, ttn )
584  use m_fstr
585  use m_matmatrix
586  use mmechgauss
587  use m_utilities
588 
589  type(tgaussstatus), intent(inout) :: gauss
590  integer(kind=kint), intent(in) :: flag
591  real(kind=kreal), intent(in) :: rot(3,3)
592  real(kind=kreal), intent(in) :: dstrain(6)
593  real(kind=kreal), intent(in) :: f(3,3) !deformation gradient (used for ss_out)
594  real(kind=kreal), intent(in) :: coordsys(3,3)
595  real(kind=kreal), intent(in) :: time
596  real(kind=kreal), intent(in) :: tincr
597  real(kind=kreal), intent(in), optional :: ttc
598  real(kind=kreal), intent(in), optional :: tt0
599  real(kind=kreal), intent(in), optional :: ttn
600 
601  integer(kind=kint) :: mtype, i, j, k
602  integer(kind=kint) :: isep
603  real(kind=kreal) :: d(6,6), dstress(6), dumstress(3,3), dum(3,3), trd, det
604  real(kind=kreal) :: tensor(6)
605  real(kind=kreal) :: eigval(3)
606  real(kind=kreal) :: princ(3,3), norm
607 
608  mtype = gauss%pMaterial%mtype
609 
610  if( iselastoplastic(mtype) .OR. mtype == norton )then
611  isep = 1
612  else
613  isep = 0
614  endif
615 
616  if( present(ttc) .AND. present(ttn) ) then
617  call matlmatrix( gauss, d3, d, time, tincr, coordsys, ttc, isep )
618  else
619  call matlmatrix( gauss, d3, d, time, tincr, coordsys, isep=isep)
620  end if
621 
622  if( flag == infinitesimal ) then
623 
624  gauss%stress(1:6) = matmul( d(1:6, 1:6), dstrain(1:6) )
625  if( isviscoelastic(mtype) .AND. tincr /= 0.0d0 ) then
626  if( present(ttc) .AND. present(ttn) ) then
627  call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr, ttc, ttn )
628  else
629  call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr )
630  end if
631  gauss%stress = real(gauss%stress)
632  end if
633 
634  else if( flag == totallag ) then
635 
636  if( ishyperelastic(mtype) .OR. mtype == userelastic .OR. mtype == usermaterial ) then
637  call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys )
638  else if( ( isviscoelastic(mtype) .OR. mtype == norton ) .AND. tincr /= 0.0d0 ) then
639  gauss%pMaterial%mtype=mtype
640  if( present(ttc) .AND. present(ttn) ) then
641  call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr, ttc, ttn )
642  else
643  call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr )
644  end if
645  else
646  gauss%stress(1:6) = matmul( d(1:6, 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  if( present(ttc) .AND. present(ttn) ) then
655  call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr, ttc, tt0 )
656  else
657  call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys, time, tincr )
658  end if
659  else
660  dstress = real( matmul( d(1:6,1:6), dstrain(1:6) ) )
661  dumstress(1,1) = gauss%stress_bak(1)
662  dumstress(2,2) = gauss%stress_bak(2)
663  dumstress(3,3) = gauss%stress_bak(3)
664  dumstress(1,2) = gauss%stress_bak(4); dumstress(2,1)=dumstress(1,2)
665  dumstress(2,3) = gauss%stress_bak(5); dumstress(3,2)=dumstress(2,3)
666  dumstress(3,1) = gauss%stress_bak(6); dumstress(1,3)=dumstress(3,1)
667 
668  !stress integration
669  trd = dstrain(1)+dstrain(2)+dstrain(3)
670  dum(:,:) = dumstress + matmul( rot,dumstress ) - matmul( dumstress, rot ) + dumstress*trd
671  !call Hughes_Winget_rotation_3D( rot, dumstress, dum )
672 
673  gauss%stress(1) = dum(1,1) + dstress(1)
674  gauss%stress(2) = dum(2,2) + dstress(2)
675  gauss%stress(3) = dum(3,3) + dstress(3)
676  gauss%stress(4) = dum(1,2) + dstress(4)
677  gauss%stress(5) = dum(2,3) + dstress(5)
678  gauss%stress(6) = dum(3,1) + dstress(6)
679 
680  if( mtype == usermaterial ) then
681  call stressupdate( gauss, d3, dstrain, gauss%stress, coordsys )
682  else if( mtype == norton ) then
683  !gauss%pMaterial%mtype = mtype
684  if( tincr /= 0.0d0 .AND. any( gauss%stress /= 0.0d0 ) ) then
685  !gauss%pMaterial%mtype = mtype
686  if( present(ttc) .AND. present(ttn) ) then
687  call stressupdate( gauss, d3, gauss%strain, gauss%stress, coordsys, time, tincr, ttc, ttn )
688  else
689  call stressupdate( gauss, d3, gauss%strain, gauss%stress, coordsys, time, tincr )
690  end if
691  end if
692  end if
693  end if
694 
695  end if
696 
697  if( iselastoplastic(mtype) ) then
698  if( present(ttc) ) then
699  call backwardeuler( gauss%pMaterial, gauss%stress, gauss%plstrain, &
700  gauss%istatus(1), gauss%fstatus, ttc )
701  else
702  call backwardeuler( gauss%pMaterial, gauss%stress, gauss%plstrain, &
703  gauss%istatus(1), gauss%fstatus )
704  end if
705  end if
706 
707  !convert stress/strain measure for output
708  if( opsstype == kopss_solution ) then
709 
710  if( flag == infinitesimal ) then !linear
711  gauss%stress_out(1:6) = gauss%stress(1:6)
712  gauss%strain_out(1:6) = gauss%strain(1:6)
713  else !nonlinear
714  !convert stress
715  if( flag == totallag ) then
716  dumstress(1,1) = gauss%stress(1)
717  dumstress(2,2) = gauss%stress(2)
718  dumstress(3,3) = gauss%stress(3)
719  dumstress(1,2) = gauss%stress(4); dumstress(2,1)=dumstress(1,2)
720  dumstress(2,3) = gauss%stress(5); dumstress(3,2)=dumstress(2,3)
721  dumstress(3,1) = gauss%stress(6); dumstress(1,3)=dumstress(3,1)
722 
723  det = determinant33(f)
724  if( det == 0.d0 ) stop "Fail to convert stress: detF=0"
725  ! cauchy stress = (1/detF)*F*(2ndPK stress)*F^T
726  dumstress(1:3,1:3) = matmul(dumstress(1:3,1:3),transpose(f(1:3,1:3)))
727  dumstress(1:3,1:3) = (1.d0/det)*matmul(f(1:3,1:3),dumstress(1:3,1:3))
728 
729  gauss%stress_out(1) = dumstress(1,1)
730  gauss%stress_out(2) = dumstress(2,2)
731  gauss%stress_out(3) = dumstress(3,3)
732  gauss%stress_out(4) = dumstress(1,2)
733  gauss%stress_out(5) = dumstress(2,3)
734  gauss%stress_out(6) = dumstress(3,1)
735  else if( flag == updatelag ) then
736  gauss%stress_out(1:6) = gauss%stress(1:6)
737  endif
738 
739  !calc logarithmic strain
740  dum(1:3,1:3) = matmul(f(1:3,1:3),transpose(f(1:3,1:3)))
741  tensor(1) = dum(1,1)
742  tensor(2) = dum(2,2)
743  tensor(3) = dum(3,3)
744  tensor(4) = dum(1,2)
745  tensor(5) = dum(2,3)
746  tensor(6) = dum(3,1)
747  call get_principal(tensor, eigval, princ)
748 
749  do k=1,3
750  if( eigval(k) <= 0.d0 ) stop "Fail to calc log strain: stretch<0"
751  eigval(k) = 0.5d0*dlog(eigval(k)) !log(sqrt(lambda))
752  norm = dsqrt(dot_product(princ(1:3,k),princ(1:3,k)))
753  if( norm <= 0.d0 ) stop "Fail to calc log strain: stretch direction vector=0"
754  princ(1:3,k) = princ(1:3,k)/norm
755  end do
756  do i=1,3
757  do j=1,3
758  dum(i,j) = 0.d0
759  do k=1,3
760  dum(i,j) = dum(i,j) + eigval(k)*princ(i,k)*princ(j,k)
761  end do
762  end do
763  end do
764  gauss%strain_out(1) = dum(1,1)
765  gauss%strain_out(2) = dum(2,2)
766  gauss%strain_out(3) = dum(3,3)
767  gauss%strain_out(4) = 2.d0*dum(1,2)
768  gauss%strain_out(5) = 2.d0*dum(2,3)
769  gauss%strain_out(6) = 2.d0*dum(3,1)
770  endif
771 
772  else
773  gauss%stress_out(1:6) = gauss%stress(1:6)
774  gauss%strain_out(1:6) = gauss%strain(1:6)
775  end if
776 
777  end subroutine
778 
780  !---------------------------------------------------------------------*
781  subroutine update_c3 &
782  (etype,nn,ecoord, u, ddu, cdsys_id, coords, qf, &
783  gausses, iter, time, tincr, tt, t0, tn)
784  !---------------------------------------------------------------------*
785 
786  use m_fstr
787  use mmaterial
788  use mmechgauss
789  use m_matmatrix
790  use m_elastoplastic
791  use m_utilities
792 
793  integer(kind=kint), intent(in) :: etype
794  integer(kind=kint), intent(in) :: nn
795  real(kind=kreal), intent(in) :: ecoord(3, nn)
796  real(kind=kreal), intent(in) :: u(3, nn)
797  real(kind=kreal), intent(in) :: ddu(3, nn)
798  integer(kind=kint), intent(in) :: cdsys_id
799  real(kind=kreal), intent(inout) :: coords(3, 3)
800  real(kind=kreal), intent(out) :: qf(nn*3)
801  type(tgaussstatus), intent(inout) :: gausses(:)
802  integer, intent(in) :: iter
803  real(kind=kreal), intent(in) :: time
804  real(kind=kreal), intent(in) :: tincr
805  real(kind=kreal), intent(in), optional :: tt(nn)
806  real(kind=kreal), intent(in), optional :: t0(nn)
807  real(kind=kreal), intent(in), optional :: tn(nn)
808 
809  ! LCOAL VARIAVLES
810  integer(kind=kint) :: flag
811  integer(kind=kint), parameter :: ndof = 3
812  real(kind=kreal) :: b(6,ndof*nn), b1(6,ndof*nn), spfunc(nn), ina(1)
813  real(kind=kreal) :: gderiv(nn,3), gderiv1(nn,3), gdispderiv(3,3), f(3,3), det, det1, wg, ttc,tt0, ttn
814  integer(kind=kint) :: j, lx, serr
815  real(kind=kreal) :: naturalcoord(3), rot(3,3), epsth(6)
816  real(kind=kreal) :: totaldisp(3,nn), elem(3,nn), elem1(3,nn), coordsys(3,3)
817  real(kind=kreal) :: dstrain(6)
818  real(kind=kreal) :: alpo(3)
819  logical :: ierr, matlaniso
820 
821  qf(:) = 0.0d0
822  ! we suppose the same material type in the element
823  flag = gausses(1)%pMaterial%nlgeom_flag
824  elem(:,:) = ecoord(:,:)
825  totaldisp(:,:) = u(:,:)+ ddu(:,:)
826  if( flag == updatelag ) then
827  elem(:,:) = (0.5d0*ddu(:,:)+u(:,:) ) +ecoord(:,:)
828  elem1(:,:) = (ddu(:,:)+u(:,:) ) +ecoord(:,:)
829  ! elem = elem1
830  totaldisp(:,:) = ddu(:,:)
831  end if
832 
833  matlaniso = .false.
834  if( present(tt) .and. cdsys_id > 0 ) then
835  ina = tt(1)
836  call fetch_tabledata( mc_orthoexp, gausses(1)%pMaterial%dict, alpo(:), ierr, ina )
837  if( .not. ierr ) matlaniso = .true.
838  end if
839 
840  do lx = 1, numofquadpoints(etype)
841 
842  call getquadpoint( etype, lx, naturalcoord(:) )
843  call getglobalderiv(etype, nn, naturalcoord, elem, det, gderiv)
844 
845  if( cdsys_id > 0 ) then
846  call set_localcoordsys( coords, g_localcoordsys(cdsys_id), coordsys(:,:), serr )
847  if( serr == -1 ) stop "Fail to setup local coordinate"
848  if( serr == -2 ) then
849  write(*, *) "WARNING! Cannot setup local coordinate, it is modified automatically"
850  end if
851  end if
852 
853  ! ========================================================
854  ! UPDATE STRAIN and STRESS
855  ! ========================================================
856 
857  ! Thermal Strain
858  epsth = 0.0d0
859  if( present(tt) .AND. present(t0) ) then
860  call getshapefunc(etype, naturalcoord, spfunc)
861  ttc = dot_product(tt, spfunc)
862  tt0 = dot_product(t0, spfunc)
863  ttn = dot_product(tn, spfunc)
864  call cal_thermal_expansion_c3( tt0, ttc, gausses(lx)%pMaterial, coordsys, matlaniso, epsth )
865  end if
866 
867  ! Update strain
868  ! Small strain
869  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
870  dstrain(1) = gdispderiv(1, 1)
871  dstrain(2) = gdispderiv(2, 2)
872  dstrain(3) = gdispderiv(3, 3)
873  dstrain(4) = ( gdispderiv(1, 2)+gdispderiv(2, 1) )
874  dstrain(5) = ( gdispderiv(2, 3)+gdispderiv(3, 2) )
875  dstrain(6) = ( gdispderiv(3, 1)+gdispderiv(1, 3) )
876  dstrain(:) = dstrain(:)-epsth(:) ! allright?
877 
878  f(1:3,1:3) = 0.d0; f(1,1)=1.d0; f(2,2)=1.d0; f(3,3)=1.d0; !deformation gradient
879  if( flag == infinitesimal ) then
880  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(1:6)
881 
882  else if( flag == totallag ) then
883  ! Green-Lagrange strain
884  dstrain(1) = dstrain(1)+0.5d0*dot_product( gdispderiv(:, 1), gdispderiv(:, 1) )
885  dstrain(2) = dstrain(2)+0.5d0*dot_product( gdispderiv(:, 2), gdispderiv(:, 2) )
886  dstrain(3) = dstrain(3)+0.5d0*dot_product( gdispderiv(:, 3), gdispderiv(:, 3) )
887  dstrain(4) = dstrain(4)+( gdispderiv(1, 1)*gdispderiv(1, 2) &
888  +gdispderiv(2, 1)*gdispderiv(2, 2)+gdispderiv(3, 1)*gdispderiv(3, 2) )
889  dstrain(5) = dstrain(5)+( gdispderiv(1, 2)*gdispderiv(1, 3) &
890  +gdispderiv(2, 2)*gdispderiv(2, 3)+gdispderiv(3, 2)*gdispderiv(3, 3) )
891  dstrain(6) = dstrain(6)+( gdispderiv(1, 1)*gdispderiv(1, 3) &
892  +gdispderiv(2, 1)*gdispderiv(2, 3)+gdispderiv(3, 1)*gdispderiv(3, 3) )
893 
894  gausses(lx)%strain(1:6) = dstrain(1:6)+epsth(:)
895  f(1:3,1:3) = f(1:3,1:3) + gdispderiv(1:3,1:3)
896 
897  else if( flag == updatelag ) then
898  rot = 0.0d0
899  rot(1, 2)= 0.5d0*(gdispderiv(1, 2)-gdispderiv(2, 1) ); rot(2, 1) = -rot(1, 2)
900  rot(2, 3)= 0.5d0*(gdispderiv(2, 3)-gdispderiv(3, 2) ); rot(3, 2) = -rot(2, 3)
901  rot(1, 3)= 0.5d0*(gdispderiv(1, 3)-gdispderiv(3, 1) ); rot(3, 1) = -rot(1, 3)
902 
903  gausses(lx)%strain(1:6) = gausses(lx)%strain_bak(1:6)+dstrain(1:6)+epsth(:)
904 
905  call getglobalderiv(etype, nn, naturalcoord, ecoord, det1, gderiv1)
906  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) )
907 
908  end if
909 
910  ! Update stress
911  if( present(tt) .AND. present(t0) ) then
912  call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr, ttc, tt0, ttn )
913  else
914  call update_stress3d( flag, gausses(lx), rot, dstrain, f, coordsys, time, tincr )
915  end if
916 
917  ! ========================================================
918  ! calculate the internal force ( equivalent nodal force )
919  ! ========================================================
920  ! Small strain
921  b(1:6, 1:nn*ndof) = 0.0d0
922  do j=1,nn
923  b(1,3*j-2) = gderiv(j, 1)
924  b(2,3*j-1) = gderiv(j, 2)
925  b(3,3*j ) = gderiv(j, 3)
926  b(4,3*j-2) = gderiv(j, 2)
927  b(4,3*j-1) = gderiv(j, 1)
928  b(5,3*j-1) = gderiv(j, 3)
929  b(5,3*j ) = gderiv(j, 2)
930  b(6,3*j-2) = gderiv(j, 3)
931  b(6,3*j ) = gderiv(j, 1)
932  end do
933 
934  ! calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
935  if( flag == infinitesimal ) then
936 
937  else if( flag == totallag ) then
938 
939  gdispderiv(1:ndof, 1:ndof) = matmul( totaldisp(1:ndof, 1:nn), gderiv(1:nn, 1:ndof) )
940  b1(1:6, 1:nn*ndof)=0.0d0
941  do j = 1,nn
942  b1(1, 3*j-2) = gdispderiv(1, 1)*gderiv(j, 1)
943  b1(1, 3*j-1) = gdispderiv(2, 1)*gderiv(j, 1)
944  b1(1, 3*j ) = gdispderiv(3, 1)*gderiv(j, 1)
945  b1(2, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 2)
946  b1(2, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 2)
947  b1(2, 3*j ) = gdispderiv(3, 2)*gderiv(j, 2)
948  b1(3, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 3)
949  b1(3, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 3)
950  b1(3, 3*j ) = gdispderiv(3, 3)*gderiv(j, 3)
951  b1(4, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 2)
952  b1(4, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 2)
953  b1(4, 3*j ) = gdispderiv(3, 2)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 2)
954  b1(5, 3*j-2) = gdispderiv(1, 2)*gderiv(j, 3)+gdispderiv(1, 3)*gderiv(j, 2)
955  b1(5, 3*j-1) = gdispderiv(2, 2)*gderiv(j, 3)+gdispderiv(2, 3)*gderiv(j, 2)
956  b1(5, 3*j ) = gdispderiv(3, 2)*gderiv(j, 3)+gdispderiv(3, 3)*gderiv(j, 2)
957  b1(6, 3*j-2) = gdispderiv(1, 3)*gderiv(j, 1)+gdispderiv(1, 1)*gderiv(j, 3)
958  b1(6, 3*j-1) = gdispderiv(2, 3)*gderiv(j, 1)+gdispderiv(2, 1)*gderiv(j, 3)
959  b1(6, 3*j ) = gdispderiv(3, 3)*gderiv(j, 1)+gdispderiv(3, 1)*gderiv(j, 3)
960  end do
961  ! BL = BL0 + BL1
962  do j=1,nn*ndof
963  b(:,j) = b(:,j)+b1(:,j)
964  end do
965 
966  else if( flag == updatelag ) then
967 
968  call getglobalderiv(etype, nn, naturalcoord, elem1, det, gderiv)
969  b(1:6, 1:nn*ndof) = 0.0d0
970  do j = 1, nn
971  b(1, 3*j-2) = gderiv(j, 1)
972  b(2, 3*j-1) = gderiv(j, 2)
973  b(3, 3*j ) = gderiv(j, 3)
974  b(4, 3*j-2) = gderiv(j, 2)
975  b(4, 3*j-1) = gderiv(j, 1)
976  b(5, 3*j-1) = gderiv(j, 3)
977  b(5, 3*j ) = gderiv(j, 2)
978  b(6, 3*j-2) = gderiv(j, 3)
979  b(6, 3*j ) = gderiv(j, 1)
980  end do
981 
982  end if
983 
984  ! calculate the Internal Force
985  wg=getweight( etype, lx )*det
986  qf(1:nn*ndof) &
987  = qf(1:nn*ndof)+matmul( gausses(lx)%stress(1:6), b(1:6,1:nn*ndof) )*wg
988 
989  end do
990 
991  end subroutine update_c3
992  !
993  !----------------------------------------------------------------------*
994  subroutine nodalstress_c3(etype, nn, gausses, ndstrain, ndstress)
995  !----------------------------------------------------------------------*
996  !
997  ! Calculate Strain and Stress increment of solid elements
998  !
999  use mmechgauss
1000 
1001  !---------------------------------------------------------------------
1002 
1003  integer(kind=kint), intent(in) :: etype, nn
1004  type(tgaussstatus), intent(in) :: gausses(:)
1005  real(kind=kreal), intent(out) :: ndstrain(nn,6)
1006  real(kind=kreal), intent(out) :: ndstress(nn,6)
1007 
1008  !---------------------------------------------------------------------
1009 
1010  integer :: i, ic
1011  real(kind=kreal) :: temp(12)
1012 
1013  !---------------------------------------------------------------------
1014 
1015  temp(:) = 0.0d0
1016 
1017  ic = numofquadpoints(etype)
1018 
1019  do i = 1, ic
1020  temp(1:6) = temp(1:6) +gausses(i)%strain_out(1:6)
1021  temp(7:12) = temp(7:12)+gausses(i)%stress_out(1:6)
1022  end do
1023 
1024  temp(1:12) = temp(1:12)/ic
1025 
1026  forall( i=1:nn )
1027  ndstrain(i, 1:6) = temp(1:6)
1028  ndstress(i, 1:6) = temp(7:12)
1029  end forall
1030 
1031  end subroutine nodalstress_c3
1032 
1033 
1034  !----------------------------------------------------------------------*
1035  subroutine elementstress_c3(etype, gausses, strain, stress)
1036  !----------------------------------------------------------------------*
1037  !
1038  ! Calculate Strain and Stress increment of solid elements
1039  !
1040  use mmechgauss
1041 
1042  !---------------------------------------------------------------------
1043 
1044  integer(kind=kint), intent(in) :: etype
1045  type(tgaussstatus), intent(in) :: gausses(:)
1046  real(kind=kreal), intent(out) :: strain(6)
1047  real(kind=kreal), intent(out) :: stress(6)
1048 
1049  !---------------------------------------------------------------------
1050 
1051  integer :: i, ic
1052 
1053  !---------------------------------------------------------------------
1054 
1055  strain(:) = 0.0d0; stress(:) = 0.0d0
1056 
1057  ic = numofquadpoints(etype)
1058 
1059  do i = 1, ic
1060  strain(:) = strain(:)+gausses(i)%strain_out(1:6)
1061  stress(:) = stress(:)+gausses(i)%stress_out(1:6)
1062  enddo
1063 
1064  strain(:) = strain(:)/ic
1065  stress(:) = stress(:)/ic
1066 
1067  end subroutine elementstress_c3
1068 
1069 
1071  !----------------------------------------------------------------------*
1072  real(kind=kreal) function volume_c3(etype, nn, XX, YY, ZZ)
1073  !----------------------------------------------------------------------*
1074 
1075  integer(kind=kint), intent(in) :: etype, nn
1076  real(kind=kreal), intent(in) :: xx(:), yy(:), zz(:)
1077 
1078  !---------------------------------------------------------------------
1079 
1080  real(kind=kreal) :: xj(3, 3), det
1081  integer(kind=kint) :: lx
1082  real(kind=kreal) :: localcoord(3), deriv(nn, 3)
1083 
1084  !---------------------------------------------------------------------
1085 
1086  volume_c3 = 0.0d0
1087 
1088  ! LOOP FOR INTEGRATION POINTS
1089  do lx = 1, numofquadpoints(etype)
1090 
1091  call getquadpoint(etype, lx, localcoord)
1092  call getshapederiv(etype, localcoord, deriv)
1093 
1094  ! JACOBI MATRIX
1095  xj(1, 1:3)= matmul( xx(1:nn), deriv(1:nn,1:3) )
1096  xj(2, 1:3)= matmul( yy(1:nn), deriv(1:nn,1:3) )
1097  xj(3, 1:3)= matmul( zz(1:nn), deriv(1:nn,1:3) )
1098 
1099  ! DETERMINANT OF JACOBIAN
1100  det = xj(1, 1)*xj(2, 2)*xj(3, 3) &
1101  +xj(2, 1)*xj(3, 2)*xj(1, 3) &
1102  +xj(3, 1)*xj(1, 2)*xj(2, 3) &
1103  -xj(3, 1)*xj(2, 2)*xj(1, 3) &
1104  -xj(2, 1)*xj(1, 2)*xj(3, 3) &
1105  -xj(1, 1)*xj(3, 2)*xj(2, 3)
1106 
1107  volume_c3 = volume_c3+getweight(etype, lx)*det
1108 
1109  end do
1110 
1111  end function volume_c3
1112 
1113 
1114 end module m_static_lib_3d
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coodinate system.
Definition: element.f90:640
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
Definition: element.f90:126
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:482
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:734
subroutine getsubface(intype, innumber, outtype, nodes)
Find the definition of surface of the element.
Definition: element.f90:188
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate deivatives of shape fucntion in natural coordiante system.
Definition: element.f90:571
real(kind=kreal) function getweight(fetype, np)
Fetch the weight value in given gauss point.
Definition: element.f90:528
integer function numofquadpoints(fetype)
Obtains the number of quadrature points of the element.
Definition: element.f90:445
real(kind=kreal) function, dimension(3) surfacenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 3d-surface.
Definition: element.f90:863
Definition: hecmw.f90:6
This modules defines common structures for fem analysis.
type(tlocalcoordsys), dimension(:), pointer, save g_localcoordsys
subroutine set_localcoordsys(coords, coordsys, outsys, ierr)
setup of coordinate system
This module provide functions for elastoplastic calculation.
This module defined coomon data and basic structures for analysis.
Definition: m_fstr.f90:15
integer(kind=kint), parameter kopss_solution
Definition: m_fstr.f90:114
integer(kind=kint) opsstype
Definition: m_fstr.f90:116
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.f90:120
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
subroutine matlmatrix(gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp)
Calculate constituive matrix.
subroutine stressupdate(gauss, sectType, strain, stress, cdsys, time, dtime, temp, tempn)
Update strain and stress for elastic and hyperelastic materials.
This module provide common functions of Solid elements.
subroutine elementstress_c3(etype, gausses, strain, stress)
subroutine tload_c3(etype, nn, XX, YY, ZZ, TT, T0, gausses, VECT, cdsys_ID, coords)
This subroutien calculate thermal loading.
subroutine dl_c3(etype, nn, XX, YY, ZZ, RHO, LTYPE, PARAMS, VECT, NSIZE)
Distrubuted external load.
subroutine hughes_winget_rotation_3d(rot, stress_in, stress_out)
subroutine geomat_c3(stress, mat)
subroutine update_c3(etype, nn, ecoord, u, ddu, cdsys_ID, coords, qf, gausses, iter, time, tincr, TT, T0, TN)
Update strain and stress inside element.
subroutine cal_thermal_expansion_c3(tt0, ttc, material, coordsys, matlaniso, EPSTH)
subroutine stf_c3(etype, nn, ecoord, gausses, stiff, cdsys_ID, coords, time, tincr, u, temperature)
This subroutine calculate stiff matrix of general solid elements.
real(kind=kreal) function volume_c3(etype, nn, XX, YY, ZZ)
Volume of element.
subroutine update_stress3d(flag, gauss, rot, dstrain, F, coordsys, time, tincr, ttc, tt0, ttn)
subroutine nodalstress_c3(etype, nn, gausses, ndstrain, ndstress)
This module provides aux functions.
Definition: utilities.f90:6
real(kind=kreal) function determinant33(XJ)
Compute determinant for 3*3 matrix.
Definition: utilities.f90:233
subroutine get_principal(tensor, eigval, princmatrix)
Definition: utilities.f90:374
subroutine transformation(jacob, tm)
Definition: utilities.f90:339
subroutine calinverse(NN, A)
calculate inverse of matrix a
Definition: utilities.f90:258
This module summarizes all infomation of material properties.
Definition: material.f90:6
integer(kind=kint), parameter totallag
Definition: material.f90:14
character(len=dict_key_length) mc_orthoexp
Definition: material.f90:123
integer(kind=kint), parameter infinitesimal
Definition: material.f90:13
integer(kind=kint), parameter updatelag
Definition: material.f90:15
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13