FrontISTR  5.7.0
Large-scale structural analysis program with finit element method
static_LIB_2d.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_2d
7  use hecmw, only : kint, kreal
8  use elementinfo
9  implicit none
10 
11 contains
12  !***********************************************************************
13  ! 2D Element:
14  ! STF_C2 :Calculate stiff matrix of 2d elements
15  ! DL_C2 :Deal with DLOAD conditions
16  ! TLOAD_C2 :Deal with thermal expansion force
17  ! UpdateST_C2 : Update Strain and stress
18  !***********************************************************************
19  !----------------------------------------------------------------------*
20  subroutine stf_c2( ETYPE,NN,ecoord,gausses,PARAM1,stiff &
21  ,ISET, u)
22  !----------------------------------------------------------------------*
23  use mmechgauss
24  use m_matmatrix
25  integer(kind=kint), intent(in) :: etype
26  integer(kind=kint), intent(in) :: nn
27  type(tgaussstatus), intent(in) :: gausses(:)
28  real(kind=kreal), intent(in) :: ecoord(2,nn),param1
29  real(kind=kreal), intent(out) :: stiff(:,:)
30  integer(kind=kint),intent(in) :: ISET
31  real(kind=kreal), intent(in), optional :: u(:,:)
32 
33  ! LOCAL VARIABLES
34  integer(kind=kint) :: flag
35  integer(kind=kint) NDOF
36  parameter(ndof=2)
37  real(kind=kreal) d(4,4),b(4,ndof*nn),db(4,ndof*nn)
38  real(kind=kreal) h(nn),stress(4)
39  real(kind=kreal) thick,pai
40  real(kind=kreal) det,rr,wg
41  real(kind=kreal) localcoord(2)
42  real(kind=kreal) gderiv(nn,2), cdsys(3,3)
43  integer(kind=kint) J,LX
44  real(kind=kreal) gdispderiv(2,2)
45  real(kind=kreal) b1(4,nn*ndof)
46  real(kind=kreal) smat(4,4)
47  real(kind=kreal) bn(4,nn*ndof), sbn(4,nn*ndof)
48 
49  stiff =0.d0
50  flag = gausses(1)%pMaterial%nlgeom_flag
51  ! THICKNESS
52  thick=param1
53  ! FOR AX-SYM. ANALYSIS
54  if(iset==2) then
55  thick=1.d0
56  pai=4.d0*atan(1.d0)
57  endif
58  !* LOOP OVER ALL INTEGRATION POINTS
59  do lx=1,numofquadpoints( etype )
60  call matlmatrix( gausses(lx), iset, d, 1.d0, 1.d0, cdsys, 0.d0 )
61  if( .not. present(u) ) flag=infinitesimal ! enforce to infinitesimal deformation analysis
62 
63  if( flag==1 .and. iset == 2 ) then
64  write(*,'(a)') ' PROGRAM STOP : non-TL element for axisymmetric element'
65  stop
66  endif
67 
68  call getquadpoint( etype, lx, localcoord )
69  call getglobalderiv( etype, nn, localcoord, ecoord, det, gderiv )
70  !
71  if(iset==2) then
72  call getshapefunc( etype, localcoord, h(:) )
73  rr=dot_product( h(1:nn), ecoord(1,1:nn) )
74  wg=getweight( etype, lx )*det*rr*2.d0*pai
75  else
76  rr=thick
77  h(:)=0.d0
78  wg=getweight( etype, lx )*det*rr
79  end if
80  do j=1,nn
81  b(1,2*j-1)=gderiv(j,1)
82  b(2,2*j-1)=0.d0
83  b(3,2*j-1)=gderiv(j,2)
84  b(1,2*j )=0.d0
85  b(2,2*j )=gderiv(j,2)
86  b(3,2*j )=gderiv(j,1)
87  b(4,2*j-1)=h(j)/rr
88  b(4,2*j )=0.d0
89  enddo
90  !
91  ! ----- calculate the BL1 matrix ( TOTAL LAGRANGE METHOD )
92  if( flag==1 ) then
93  ! ----- dudx(i,j) ==> gdispderiv(i,j)
94  gdispderiv(1:ndof, 1:ndof) = matmul( u(1:ndof,1:nn), gderiv(1:nn,1:ndof) )
95  b1(1:4,1:nn*ndof) = 0.d0
96  do j=1,nn
97  b1(1,2*j-1) = gdispderiv(1,1)*gderiv(j,1)
98  b1(2,2*j-1) = gdispderiv(1,2)*gderiv(j,2)
99  b1(3,2*j-1) = gdispderiv(1,2)*gderiv(j,1)+gdispderiv(1,1)*gderiv(j,2)
100  b1(1,2*j ) = gdispderiv(2,1)*gderiv(j,1)
101  b1(2,2*j ) = gdispderiv(2,2)*gderiv(j,2)
102  b1(3,2*j ) = gdispderiv(2,2)*gderiv(j,1)+gdispderiv(2,1)*gderiv(j,2)
103  b1(4,2*j-1) = 0.d0
104  b1(4,2*j ) = 0.d0
105  enddo
106  do j=1,nn*ndof
107  b(:,j) = b(:,j)+b1(:,j)
108  enddo
109  endif
110  !
111  db(1:4,1:nn*ndof) = 0.d0
112  db(1:4,1:nn*ndof) = matmul( d(1:4,1:4), b(1:4,1:nn*ndof) )
113  stiff(1:nn*ndof,1:nn*ndof) = stiff(1:nn*ndof,1:nn*ndof) + &
114  matmul( transpose( b(1:4,1:nn*ndof) ), db(1:4,1:nn*ndof) )*wg
115  !
116  ! ----- calculate the stress matrix ( TOTAL LAGRANGE METHOD )
117  if( flag==1 ) then
118  stress(1:4)=gausses(lx)%stress(1:4)
119  bn(1:4,1:nn*ndof) = 0.d0
120  do j=1,nn
121  bn(1,2*j-1) = gderiv(j,1)
122  bn(2,2*j ) = gderiv(j,1)
123  bn(3,2*j-1) = gderiv(j,2)
124  bn(4,2*j ) = gderiv(j,2)
125  enddo
126  smat(:,:) = 0.d0
127  do j=1,2
128  smat(j ,j ) = stress(1)
129  smat(j ,j+2) = stress(3)
130  smat(j+2,j ) = stress(3)
131  smat(j+2,j+2) = stress(2)
132  enddo
133  sbn(1:4,1:nn*ndof) = matmul( smat(1:4,1:4), bn(1:4,1:nn*ndof) )
134  stiff(1:nn*ndof,1:nn*ndof) = stiff(1:nn*ndof,1:nn*ndof) + &
135  matmul( transpose( bn(1:4,1:nn*ndof) ), sbn(1:4,1:nn*ndof) )*wg
136  endif
137  !
138  enddo
139 
140  end subroutine stf_c2
141 
142 
143  !----------------------------------------------------------------------*
144  subroutine dl_c2(ETYPE,NN,XX,YY,RHO,PARAM1,LTYPE,PARAMS,VECT,NSIZE,ISET)
145  !----------------------------------------------------------------------*
146  !**
147  !** Deal with DLOAD conditions
148  !**
149  ! BX LTYPE=1 :BODY FORCE IN X-DIRECTION
150  ! BY LTYPE=2 :BODY FORCE IN Y-DIRECTION
151  ! GRAV LTYPE=4 :GRAVITY FORCE
152  ! CENT LTYPE=5 :CENTRIFUGAL LOAD
153  ! P1 LTYPE=10 :TRACTION IN NORMAL-DIRECTION FOR FACE-1
154  ! P2 LTYPE=20 :TRACTION IN NORMAL-DIRECTION FOR FACE-2
155  ! P3 LTYPE=30 :TRACTION IN NORMAL-DIRECTION FOR FACE-3
156  ! P4 LTYPE=40 :TRACTION IN NORMAL-DIRECTION FOR FACE-4
157  ! I/F VARIABLES
158  integer(kind=kint), intent(in) :: ETYPE,NN
159  real(kind=kreal), intent(in) :: xx(nn),yy(nn),params(0:6)
160  real(kind=kreal), intent(out) :: vect(nn*2)
161  real(kind=kreal) rho,param1
162  integer(kind=kint) LTYPE,NSIZE,ISET
163  ! LOCAL VARIABLES
164  integer(kind=kint), parameter :: NDOF=2
165  real(kind=kreal) h(nn)
166  real(kind=kreal) plx(nn),ply(nn)
167  real(kind=kreal) xj(2,2),det,rr,wg
168  real(kind=kreal) pai,val,thick
169  integer(kind=kint) IVOL,ISUF
170  real(kind=kreal) ax,ay,rx,ry
171  real(kind=kreal) normal(2)
172  real(kind=kreal) coefx,coefy,xcod,ycod,hx,hy,phx,phy
173  integer(kind=kint) NOD(NN)
174  integer(kind=kint) LX,I,SURTYPE,NSUR
175  real(kind=kreal) vx,vy,localcoord(2),deriv(nn,2),elecoord(2,nn)
176 
177  rx = 0.0d0; ry = 0.0d0
178  ay = 0.0d0; ax = 0.0d0
179 
180  ! CONSTANT
181  pai=4.d0*atan(1.d0)
182  ! SET VALUE
183  val=params(0)
184  ! THICKNESS
185  thick=param1
186  ! CLEAR VECT
187  nsize=nn*ndof
188  vect(1:nsize)=0.d0
189  !
190  ! SELECTION OF LOAD TYPE
191  !
192  ivol=0
193  isuf=0
194  if( ltype.LT.10 ) then
195  ivol=1
196  if( ltype.EQ.5 ) then
197  ax=params(1)
198  ay=params(2)
199  rx=params(4)
200  ry=params(5)
201  endif
202  else if( ltype.GE.10 ) then
203  isuf=1
204  call getsubface( etype, ltype/10, surtype, nod )
205  nsur = getnumberofnodes( surtype )
206  endif
207  !** SURFACE LOAD
208  if( isuf==1 ) then
209  do i=1,nsur
210  elecoord(1,i)=xx(nod(i))
211  elecoord(2,i)=yy(nod(i))
212  enddo
213  !** INTEGRATION OVER SURFACE
214  do lx=1,numofquadpoints( surtype )
215  call getquadpoint( surtype, lx, localcoord(1:1) )
216  call getshapefunc( surtype, localcoord(1:1), h(1:nsur) )
217  normal=edgenormal( surtype, nsur, localcoord(1:1), elecoord(:,1:nsur) )
218  wg = getweight( surtype, lx )
219  if( iset==2 ) then
220  rr=0.d0
221  do i=1,nsur
222  rr=rr+h(i)*xx(nod(i))
223  enddo
224  wg=wg*rr*2.d0*pai
225  else
226  wg=wg*thick
227  endif
228  do i=1,nsur
229  vect(2*nod(i)-1)=vect(2*nod(i)-1)+val*wg*h(i)*normal(1)
230  vect(2*nod(i) )=vect(2*nod(i) )+val*wg*h(i)*normal(2)
231  enddo
232  enddo
233  endif
234  !** VOLUME LOAD
235  if( ivol==1 ) then
236  plx(:)=0.d0
237  ply(:)=0.d0
238  do lx=1,numofquadpoints( etype )
239  call getquadpoint( etype, lx, localcoord(1:2) )
240  call getshapederiv( etype, localcoord(1:2), deriv )
241  call getshapefunc( etype, localcoord(1:2), h(1:nn) )
242  xj(1,1:2)=matmul( xx(1:nn), deriv(1:nn,1:2) )
243  xj(2,1:2)=matmul( yy(1:nn), deriv(1:nn,1:2) )
244 
245  det=xj(1,1)*xj(2,2)-xj(2,1)*xj(1,2)
246 
247  wg = getweight( etype, lx )
248  if(iset==2) then
249  rr=dot_product( h(1:nn),xx(1:nn) )
250  wg=wg*det*rr*2.d0*pai
251  else
252  rr=thick
253  wg=wg*det*rr
254  end if
255  coefx=1.d0
256  coefy=1.d0
257  ! CENTRIFUGAL LOAD
258  if( ltype==5 ) then
259  xcod=dot_product( h(1:nn),xx(1:nn) )
260  ycod=dot_product( h(1:nn),yy(1:nn) )
261  hx=ax + ( (xcod-ax)*rx+(ycod-ay)*ry )/(rx**2+ry**2)*rx
262  hy=ay + ( (xcod-ax)*rx+(ycod-ay)*ry )/(rx**2+ry**2)*ry
263  phx=xcod-hx
264  phy=ycod-hy
265  coefx=rho*val*val*phx
266  coefy=rho*val*val*phy
267  end if
268  do i=1,nn
269  plx(i)=plx(i)+h(i)*wg*coefx
270  ply(i)=ply(i)+h(i)*wg*coefy
271  enddo
272  enddo
273  if( ltype.EQ.1) then
274  do i=1,nn
275  vect(2*i-1)=val*plx(i)
276  vect(2*i )=0.d0
277  enddo
278  else if( ltype.EQ.2 ) then
279  do i=1,nn
280  vect(2*i-1)=0.d0
281  vect(2*i )=val*ply(i)
282  enddo
283  else if( ltype.EQ.4 ) then
284  vx=params(1)
285  vy=params(2)
286  vx=vx/sqrt( params(1)**2 + params(2)**2 )
287  vy=vy/sqrt( params(1)**2 + params(2)**2 )
288  do i=1,nn
289  vect(2*i-1)=val*plx(i)*rho*vx
290  vect(2*i )=val*ply(i)*rho*vy
291  enddo
292  else if( ltype==5 ) then
293  do i=1,nn
294  vect(2*i-1)=plx(i)
295  vect(2*i )=ply(i)
296  enddo
297  end if
298  endif
299  end subroutine dl_c2
300 
301  !----------------------------------------------------------------------*
302  subroutine tload_c2(ETYPE,NN,XX,YY,TT,T0,gausses,PARAM1,ISET,VECT)
303  !----------------------------------------------------------------------*
304  !
305  ! THERMAL LOAD OF PLANE ELEMENT
306  !
307  use mmaterial
308  use m_matmatrix
309  use m_elasticlinear
310  ! I/F VARIABLES
311  integer(kind=kint), intent(in) :: ETYPE,NN
312  real(kind=kreal),intent(in) :: xx(nn),yy(nn),tt(nn),t0(nn),param1
313  type(tgaussstatus), intent(in) :: gausses(:)
314  real(kind=kreal),intent(out) :: vect(nn*2)
315  integer(kind=kint) iset
316  ! LOCAL VARIABLES
317  integer(kind=kint) ndof
318  parameter(ndof=2)
319  real(kind=kreal) alp,pp,d(4,4),b(4,ndof*nn)
320  real(kind=kreal) h(nn)
321  real(kind=kreal) eps(4),sgm(4),localcoord(2)
322  real(kind=kreal) deriv(nn,2),dnde(2,nn)
323  real(kind=kreal) xj(2,2),det,rr,dum
324  real(kind=kreal) xji(2,2)
325  real(kind=kreal) pai,thick,wg
326  integer(kind=kint) j,lx
327  real(kind=kreal) tempc,temp0,thermal_eps
328  !*************************
329  ! CONSTANT
330  !*************************
331  pai=4.d0*atan(1.d0)
332  ! CLEAR VECT
333  vect(1:nn*ndof)=0.d0
334  ! THICKNESS
335  thick=param1
336  ! FOR AX-SYM. ANALYSIS
337  if(iset==2) thick=1.d0
338  ! We suppose material properties doesn't varies inside element
339  call calelasticmatrix( gausses(1)%pMaterial, iset, d, 0.d0 )
340  alp = gausses(1)%pMaterial%variables(m_exapnsion)
341  pp = gausses(1)%pMaterial%variables(m_poisson)
342  !* LOOP OVER ALL INTEGRATION POINTS
343  do lx=1,numofquadpoints( etype )
344  call getquadpoint( etype, lx, localcoord )
345  call getshapefunc( etype, localcoord, h(1:nn) )
346  call getshapederiv( etype, localcoord, deriv(:,:) )
347  xj(1,1:2)=matmul( xx(1:nn), deriv(1:nn,1:2) )
348  xj(2,1:2)=matmul( yy(1:nn), deriv(1:nn,1:2) )
349 
350  det=xj(1,1)*xj(2,2)-xj(2,1)*xj(1,2)
351 
352  tempc=dot_product(h(1:nn),tt(1:nn))
353  temp0=dot_product(h(1:nn),t0(1:nn))
354  if(iset==2) then
355  rr=dot_product(h(1:nn),xx(1:nn))
356  wg=getweight( etype, lx )*det*rr*2.d0*pai
357  else
358  rr=thick
359  wg=getweight( etype, lx )*det*rr
360  end if
361  dum=1.d0/det
362  xji(1,1)= xj(2,2)*dum
363  xji(1,2)=-xj(2,1)*dum
364  xji(2,1)=-xj(1,2)*dum
365  xji(2,2)= xj(1,1)*dum
366 
367  dnde=matmul( xji, transpose(deriv) )
368  do j=1,nn
369  b(1,2*j-1)=dnde(1,j)
370  b(2,2*j-1)=0.d0
371  b(3,2*j-1)=dnde(2,j)
372  b(1,2*j )=0.d0
373  b(2,2*j )=dnde(2,j)
374  b(3,2*j )=dnde(1,j)
375  b(4,2*j-1)=0.d0
376  b(4,2*j )=0.d0
377  if( iset==2 ) then
378  b(4,2*j-1)=h(j)/rr
379  endif
380  enddo
381  !**
382  !** THERMAL EPS
383  !**
384  thermal_eps=alp*(tempc-temp0)
385  if( iset .EQ. 2 ) then
386  eps(1)=thermal_eps
387  eps(2)=thermal_eps
388  eps(3)=0.d0
389  eps(4)=thermal_eps
390  else if( iset.EQ.0 ) then
391  eps(1)=thermal_eps*(1.d0+pp)
392  eps(2)=thermal_eps*(1.d0+pp)
393  eps(3)=0.d0
394  eps(4)=0.d0
395  else if( iset.EQ.1 ) then
396  eps(1)=thermal_eps
397  eps(2)=thermal_eps
398  eps(3)=0.d0
399  eps(4)=0.d0
400  endif
401  !**
402  !** SET SGM {s}=[D]{e}
403  !**
404  sgm = matmul( eps(1:4), d(1:4,1:4) )
405  !**
406  !** CALCULATE LOAD {F}=[B]T{e}
407  !**
408  vect(1:nn*ndof)=vect(1:nn*ndof)+matmul( sgm(1:4), b(1:4,1:nn*ndof) )*wg
409  enddo
410 
411  end subroutine tload_c2
412  !
413  !
415  !---------------------------------------------------------------------*
416  subroutine update_c2( etype,nn,ecoord,gausses,PARAM1,iset, &
417  u, ddu, qf, TT, T0, TN )
418  !---------------------------------------------------------------------*
419  use m_fstr
420  use mmechgauss
421  use m_matmatrix
422  ! I/F VARIABLES
423  integer(kind=kint), intent(in) :: etype, nn
424  real(kind=kreal), intent(in) :: ecoord(3,nn),param1
425  integer(kind=kint), intent(in) :: iset
426  type(tgaussstatus), intent(inout) :: gausses(:)
427  real(kind=kreal), intent(in) :: u(2,nn)
428  real(kind=kreal), intent(in) :: ddu(2,nn)
429  real(kind=kreal), intent(out) :: qf(:)
430  real(kind=kreal), intent(in) :: tt(nn)
431  real(kind=kreal), intent(in) :: t0(nn)
432  real(kind=kreal), intent(in) :: tn(nn)
433 
434 
435  integer(kind=kint), parameter :: ndof=2
436  real(kind=kreal) :: d(4,4), b(4,ndof*nn)
437  real(kind=kreal) :: h(nn)
438  real(kind=kreal) :: thick,pai,rr
439  real(kind=kreal) :: det, wg
440  real(kind=kreal) :: localcoord(2)
441  real(kind=kreal) :: gderiv(nn,2), ttc, tt0, ttn
442  real(kind=kreal) :: cdsys(3,3)
443  integer(kind=kint) :: j, LX
444  real(kind=kreal) :: totaldisp(2*nn)
445  real(kind=kreal) :: epsth(4), alp, alp0, ina(1), outa(1)
446  logical :: ierr
447  !
448  qf(:) = 0.d0
449  do j=1,nn
450  totaldisp((j-1)*2+1) = u(1,j) + ddu(1,j)
451  totaldisp(j*2) = u(2,j) + ddu(2,j)
452  enddo
453  !
454  ! THICKNESS
455  thick=param1
456  ! FOR AX-SYM. ANALYSIS
457  if(iset==2) then
458  thick=1.d0
459  pai=4.d0*atan(1.d0)
460  endif
461  !* LOOP OVER ALL INTEGRATION POINTS
462  do lx=1, numofquadpoints(etype)
463  call matlmatrix( gausses(lx), iset, d, 1.d0, 1.d0, cdsys, 0.d0 )
464 
465  call getquadpoint( etype, lx, localcoord(:) )
466  call getglobalderiv( etype, nn, localcoord, ecoord, det, gderiv )
467  !
468  epsth = 0.d0
469  call getshapefunc( etype, localcoord, h(:) )
470  ttc = dot_product(tt(:), h(:))
471  tt0 = dot_product(t0(:), h(:))
472  ttn = dot_product(tn(:), h(:))
473  if( dabs(ttc-tt0) > 1.d-14 ) then
474  ina(1) = ttc
475  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
476  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
477  alp = outa(1)
478  ina(1) = tt0
479  call fetch_tabledata( mc_themoexp, gausses(lx)%pMaterial%dict, outa(:), ierr, ina )
480  if( ierr ) outa(1) = gausses(lx)%pMaterial%variables(m_exapnsion)
481  alp0 = outa(1)
482  epsth(1:2)=alp*(ttc-ref_temp)-alp0*(tt0-ref_temp)
483  end if
484  !
485  wg=getweight( etype, lx )*det
486  if(iset==2) then
487  call getshapefunc( etype, localcoord, h(:) )
488  rr=dot_product( h(1:nn), ecoord(1,1:nn) )
489  else
490  rr=thick
491  h(:)=0.d0
492  end if
493  do j=1,nn
494  b(1,2*j-1)=gderiv(j,1)
495  b(2,2*j-1)=0.d0
496  b(3,2*j-1)=gderiv(j,2)
497  b(1,2*j )=0.d0
498  b(2,2*j )=gderiv(j,2)
499  b(3,2*j )=gderiv(j,1)
500  b(4,2*j-1)=h(j)/rr
501  b(4,2*j )=0.d0
502  enddo
503 
504  gausses(lx)%strain(1:4) = matmul( b(:,:), totaldisp )
505  gausses(lx)%stress(1:4) = matmul( d(1:4, 1:4), gausses(lx)%strain(1:4)-epsth(1:4) )
506 
507  !set stress and strain for output
508  gausses(lx)%strain_out(1:4) = gausses(lx)%strain(1:4)
509  gausses(lx)%stress_out(1:4) = gausses(lx)%stress(1:4)
510 
511  !
512  ! ----- calculate the Internal Force
513  qf(1:nn*ndof) = qf(1:nn*ndof) + &
514  matmul( gausses(lx)%stress(1:4), b(1:4,1:nn*ndof) )*wg
515  !
516  ! calculate strain energy
517  gausses(lx)%strain_energy = 0.5d0*dot_product(gausses(lx)%stress(1:6),gausses(lx)%strain(1:6))*wg
518  enddo
519  !
520  end subroutine update_c2
521  !
522  !----------------------------------------------------------------------*
523  subroutine nodalstress_c2(ETYPE,NN,gausses,ndstrain,ndstress)
524  !----------------------------------------------------------------------*
525  !
526  ! Calculate Strain and Stress increment of solid elements
527  !
528  use mmechgauss
529 
530  integer(kind=kint), intent(in) :: ETYPE,NN
531  type(tgaussstatus), intent(in) :: gausses(:)
532  real(kind=kreal), intent(out) :: ndstrain(nn,4)
533  real(kind=kreal), intent(out) :: ndstress(nn,4)
534 
535  integer :: i,ic
536  real(kind=kreal) :: temp(8)
537 
538  temp(:)=0.d0
539  ic = numofquadpoints(etype)
540  do i=1,ic
541  temp(1:4) = temp(1:4) + gausses(i)%strain_out(1:4)
542  temp(5:8) = temp(5:8) + gausses(i)%stress_out(1:4)
543  enddo
544  temp(1:8) = temp(1:8)/ic
545  do i=1,nn
546  ndstrain(i,1:4) = temp(1:4)
547  ndstress(i,1:4) = temp(5:8)
548  enddo
549 
550  end subroutine
551 
552  !----------------------------------------------------------------------*
553  subroutine elementstress_c2(ETYPE,gausses,strain,stress)
554  !----------------------------------------------------------------------*
555  !
556  ! Calculate Strain and Stress increment of solid elements
557  !
558  use mmechgauss
559  integer(kind=kint), intent(in) :: ETYPE
560  type(tgaussstatus), intent(in) :: gausses(:)
561  real(kind=kreal), intent(out) :: strain(4)
562  real(kind=kreal), intent(out) :: stress(4)
563 
564  integer :: i,ic
565 
566  strain(:)=0.d0; stress(:)=0.d0
567  ic = numofquadpoints(etype)
568  do i=1,ic
569  strain(:) = strain(:) + gausses(i)%strain_out(1:4)
570  stress(:) = stress(:) + gausses(i)%stress_out(1:4)
571  enddo
572  strain(:) = strain(:)/ic
573  stress(:) = stress(:)/ic
574 
575  end subroutine
576 
577 end module m_static_lib_2d
elementinfo::getshapederiv
subroutine getshapederiv(fetype, localcoord, shapederiv)
Calculate derivatives of shape function in natural coordinate system.
Definition: element.f90:578
m_static_lib_2d::update_c2
subroutine update_c2(etype, nn, ecoord, gausses, PARAM1, iset, u, ddu, qf, TT, T0, TN)
Update strain and stress inside element.
Definition: static_LIB_2d.f90:418
mmaterial::m_exapnsion
integer(kind=kint), parameter m_exapnsion
Definition: material.f90:105
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
m_static_lib_2d::stf_c2
subroutine stf_c2(ETYPE, NN, ecoord, gausses, PARAM1, stiff, ISET, u)
Definition: static_LIB_2d.f90:22
mmaterial::m_poisson
integer(kind=kint), parameter m_poisson
Definition: material.f90:93
m_static_lib_2d
This module provides common functions of Plane deformation elements.
Definition: static_LIB_2d.f90:6
m_fstr::eps
real(kind=kreal) eps
Definition: m_fstr.f90:142
mmechgauss
This modules defines a structure to record history dependent parameter in static analysis.
Definition: mechgauss.f90:6
m_static_lib_2d::nodalstress_c2
subroutine nodalstress_c2(ETYPE, NN, gausses, ndstrain, ndstress)
Definition: static_LIB_2d.f90:524
elementinfo::edgenormal
real(kind=kreal) function, dimension(2) edgenormal(fetype, nn, localcoord, elecoord)
Calculate normal of 2d-edge.
Definition: element.f90:905
elementinfo::getglobalderiv
subroutine getglobalderiv(fetype, nn, localcoord, elecoord, det, gderiv)
Calculate shape derivative in global coordinate system.
Definition: element.f90:741
m_fstr
This module defines common data and basic structures for analysis.
Definition: m_fstr.f90:15
m_static_lib_2d::tload_c2
subroutine tload_c2(ETYPE, NN, XX, YY, TT, T0, gausses, PARAM1, ISET, VECT)
Definition: static_LIB_2d.f90:303
m_fstr::ref_temp
real(kind=kreal), pointer ref_temp
REFTEMP.
Definition: m_fstr.f90:136
m_elasticlinear
This module provides functions for elastic material.
Definition: ElasticLinear.f90:6
elementinfo
This module encapsulate the basic functions of all elements provide by this software.
Definition: element.f90:43
elementinfo::getquadpoint
subroutine getquadpoint(fetype, np, pos)
Fetch the coordinate of gauss point.
Definition: element.f90:489
hecmw
Definition: hecmw.f90:6
m_static_lib_2d::elementstress_c2
subroutine elementstress_c2(ETYPE, gausses, strain, stress)
Definition: static_LIB_2d.f90:554
elementinfo::getshapefunc
subroutine getshapefunc(fetype, localcoord, func)
Calculate the shape function in natural coordinate system.
Definition: element.f90:647
mmechgauss::tgaussstatus
All data should be recorded in every quadrature points.
Definition: mechgauss.f90:13
elementinfo::getsubface
subroutine getsubface(intype, innumber, outtype, nodes)
Find the definition of surface of the element.
Definition: element.f90:193
m_static_lib_2d::dl_c2
subroutine dl_c2(ETYPE, NN, XX, YY, RHO, PARAM1, LTYPE, PARAMS, VECT, NSIZE, ISET)
Definition: static_LIB_2d.f90:145
m_matmatrix
This module manages calculation relates with materials.
Definition: calMatMatrix.f90:6
mmaterial
This module summarizes all information of material properties.
Definition: material.f90:6
m_elasticlinear::calelasticmatrix
subroutine calelasticmatrix(matl, sectType, D, temp, hdflag)
Calculate isotropic elastic matrix.
Definition: ElasticLinear.f90:16
m_matmatrix::matlmatrix
subroutine matlmatrix(gauss, sectType, matrix, time, dtime, cdsys, temperature, isEp, hdflag)
Calculate constituive matrix.
Definition: calMatMatrix.f90:30
elementinfo::getnumberofnodes
integer(kind=kind(2)) function getnumberofnodes(etype)
Obtain number of nodes of the element.
Definition: element.f90:131